{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## TAnoGAN: Time Series Anomaly Detection with Generative Adversarial Networks\n",
    "\n",
    "Published in: 2020 IEEE Symposium Series on Computational Intelligence (SSCI)\n",
    "\n",
    "Paper link: https://ieeexplore.ieee.org/abstract/document/9308512"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import required libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import random\n",
    "import torch\n",
    "import torch.nn as nn\n",
    "import torch.backends.cudnn as cudnn\n",
    "import torch.optim as optim\n",
    "import torch.utils.data\n",
    "import torchvision\n",
    "import torch.nn.init as init\n",
    "from torch.autograd import Variable\n",
    "import datetime\n",
    "from nab_dataset import NabDataset\n",
    "from models.recurrent_models_pyramid import LSTMGenerator, LSTMDiscriminator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Define Basic Settings for Adversarial Training"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "class ArgsTrn:\n",
    "    workers=4\n",
    "    batch_size=32\n",
    "    epochs=20\n",
    "    lr=0.0002\n",
    "    cuda = True\n",
    "    manualSeed=2\n",
    "    \n",
    "opt_trn=ArgsTrn()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "torch.manual_seed(opt_trn.manualSeed)\n",
    "cudnn.benchmark = True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setup Data Loader"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# location of datasets and category\n",
    "end_name = 'ambient_temperature_system_failure.csv' # dataset name\n",
    "data_file = 'data\\\\realKnownCause\\\\'+end_name # dataset category and dataset name\n",
    "key = 'realKnownCause/'+end_name # This key is used for reading anomaly labels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# settings for data loader\n",
    "class DataSettings:\n",
    "    \n",
    "    def __init__(self):\n",
    "        self.BASE = 'D:\\\\ResearchDataGtx1060\\\\AnomalyDetectionData\\\\NabDataset\\\\'\n",
    "        self.label_file = 'labels\\\\combined_windows.json'\n",
    "        self.data_file = data_file\n",
    "        self.key = key\n",
    "        self.train = True\n",
    "        self.window_length = 60\n",
    "    \n",
    "    \n",
    "data_settings = DataSettings()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(torch.Size([7207, 60, 1]), torch.Size([7207]))"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# define dataset object and data loader object for NAB dataset\n",
    "dataset = NabDataset(data_settings=data_settings)\n",
    "dataloader = torch.utils.data.DataLoader(dataset, batch_size=opt_trn.batch_size,\n",
    "                                         shuffle=True, num_workers=int(opt_trn.workers))\n",
    "\n",
    "dataset.x.shape, dataset.y.shape # check the dataset shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setup Models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "device = torch.device(\"cuda:0\" if opt_trn.cuda else \"cpu\") # select the device\n",
    "seq_len = dataset.window_length # sequence length is equal to the window length\n",
    "in_dim = dataset.n_feature # input dimension is same as number of feature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create generator and discriminator models\n",
    "netD = LSTMDiscriminator(in_dim=in_dim, device=device).to(device)\n",
    "netG = LSTMGenerator(in_dim=in_dim, out_dim=in_dim, device=device).to(device)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "|Discriminator Architecture|\n",
      " LSTMDiscriminator(\n",
      "  (lstm): LSTM(1, 100, batch_first=True)\n",
      "  (linear): Sequential(\n",
      "    (0): Linear(in_features=100, out_features=1, bias=True)\n",
      "    (1): Sigmoid()\n",
      "  )\n",
      ")\n",
      "|Generator Architecture|\n",
      " LSTMGenerator(\n",
      "  (lstm0): LSTM(1, 32, batch_first=True)\n",
      "  (lstm1): LSTM(32, 64, batch_first=True)\n",
      "  (lstm2): LSTM(64, 128, batch_first=True)\n",
      "  (linear): Sequential(\n",
      "    (0): Linear(in_features=128, out_features=1, bias=True)\n",
      "    (1): Tanh()\n",
      "  )\n",
      ")\n"
     ]
    }
   ],
   "source": [
    "print(\"|Discriminator Architecture|\\n\", netD)\n",
    "print(\"|Generator Architecture|\\n\", netG)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Setup loss function\n",
    "criterion = nn.BCELoss().to(device)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# setup optimizer\n",
    "optimizerD = optim.Adam(netD.parameters(), lr=opt_trn.lr)\n",
    "optimizerG = optim.Adam(netG.parameters(), lr=opt_trn.lr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adversarial Training of Generator and Discriminator Models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-12-1b87bde9bad8>:24: UserWarning: nn.init.normal is now deprecated in favor of nn.init.normal_.\n",
      "  noise = Variable(init.normal(torch.Tensor(batch_size,seq_len,in_dim),mean=0,std=0.1)).cuda()\n",
      "<ipython-input-12-1b87bde9bad8>:38: UserWarning: nn.init.normal is now deprecated in favor of nn.init.normal_.\n",
      "  noise = Variable(init.normal(torch.Tensor(batch_size,seq_len,in_dim),mean=0,std=0.1)).cuda()\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0/20][225/226] Loss_D: 0.8151 Loss_G: 0.7403 D(x): 0.8784 D(G(z)): 0.4792 / 0.4771\n",
      "[1/20][225/226] Loss_D: 0.5119 Loss_G: 1.1240 D(x): 0.9274 D(G(z)): 0.3324 / 0.3321\n",
      "[2/20][225/226] Loss_D: 0.9016 Loss_G: 1.1996 D(x): 0.6857 D(G(z)): 0.3100 / 0.3101\n",
      "[3/20][225/226] Loss_D: 0.6742 Loss_G: 1.1108 D(x): 0.8341 D(G(z)): 0.3335 / 0.3333\n",
      "[4/20][225/226] Loss_D: 0.3134 Loss_G: 1.6948 D(x): 0.9487 D(G(z)): 0.2066 / 0.2059\n",
      "[5/20][225/226] Loss_D: 0.4672 Loss_G: 1.3923 D(x): 0.8931 D(G(z)): 0.2560 / 0.2561\n",
      "[6/20][225/226] Loss_D: 0.3287 Loss_G: 1.6844 D(x): 0.9506 D(G(z)): 0.2159 / 0.2156\n",
      "[7/20][225/226] Loss_D: 0.2732 Loss_G: 2.1543 D(x): 0.9549 D(G(z)): 0.1694 / 0.1694\n",
      "[8/20][225/226] Loss_D: 0.2306 Loss_G: 2.3086 D(x): 0.9518 D(G(z)): 0.1386 / 0.1389\n",
      "[9/20][225/226] Loss_D: 0.2907 Loss_G: 1.6242 D(x): 0.9666 D(G(z)): 0.2119 / 0.2117\n",
      "[10/20][225/226] Loss_D: 1.3005 Loss_G: 0.5467 D(x): 0.6754 D(G(z)): 0.5861 / 0.5854\n",
      "[11/20][225/226] Loss_D: 1.2143 Loss_G: 0.5787 D(x): 0.7236 D(G(z)): 0.5713 / 0.5707\n",
      "[12/20][225/226] Loss_D: 0.4149 Loss_G: 1.4677 D(x): 0.8941 D(G(z)): 0.2413 / 0.2414\n",
      "[13/20][225/226] Loss_D: 1.3836 Loss_G: 0.3552 D(x): 0.8893 D(G(z)): 0.7075 / 0.7052\n",
      "[14/20][225/226] Loss_D: 0.6905 Loss_G: 0.9112 D(x): 0.8842 D(G(z)): 0.4043 / 0.4040\n",
      "[15/20][225/226] Loss_D: 0.7800 Loss_G: 0.9786 D(x): 0.8054 D(G(z)): 0.3790 / 0.3786\n",
      "[16/20][225/226] Loss_D: 0.6193 Loss_G: 1.1315 D(x): 0.8711 D(G(z)): 0.3309 / 0.3308\n",
      "[17/20][225/226] Loss_D: 0.5345 Loss_G: 1.3347 D(x): 0.8585 D(G(z)): 0.2798 / 0.2868\n",
      "[18/20][225/226] Loss_D: 1.0625 Loss_G: 0.9381 D(x): 0.6036 D(G(z)): 0.3907 / 0.3944\n",
      "[19/20][225/226] Loss_D: 1.2283 Loss_G: 0.6555 D(x): 0.6927 D(G(z)): 0.5291 / 0.5289\n"
     ]
    }
   ],
   "source": [
    "real_label = 1\n",
    "fake_label = 0\n",
    "\n",
    "for epoch in range(opt_trn.epochs):\n",
    "    for i, (x,y) in enumerate(dataloader, 0):\n",
    "        \n",
    "        ############################\n",
    "        # (1) Update D network: maximize log(D(x)) + log(1 - D(G(z)))\n",
    "        ###########################\n",
    "\n",
    "        #Train with real data\n",
    "        netD.zero_grad()\n",
    "        real = x.to(device)\n",
    "        batch_size, seq_len = real.size(0), real.size(1)\n",
    "        label = torch.full((batch_size, seq_len, 1), real_label, device=device)\n",
    "\n",
    "        output,_ = netD.forward(real)\n",
    "        errD_real = criterion(output, label.float())\n",
    "        errD_real.backward()\n",
    "        optimizerD.step()\n",
    "        D_x = output.mean().item()\n",
    "        \n",
    "        #Train with fake data\n",
    "        noise = Variable(init.normal(torch.Tensor(batch_size,seq_len,in_dim),mean=0,std=0.1)).cuda()\n",
    "        fake,_ = netG.forward(noise)\n",
    "        output,_ = netD.forward(fake.detach()) # detach causes gradient is no longer being computed or stored to save memeory\n",
    "        label.fill_(fake_label)\n",
    "        errD_fake = criterion(output, label.float())\n",
    "        errD_fake.backward()\n",
    "        D_G_z1 = output.mean().item()\n",
    "        errD = errD_real + errD_fake\n",
    "        optimizerD.step()\n",
    "        \n",
    "        ############################\n",
    "        # (2) Update G network: maximize log(D(G(z)))\n",
    "        ###########################\n",
    "        netG.zero_grad()\n",
    "        noise = Variable(init.normal(torch.Tensor(batch_size,seq_len,in_dim),mean=0,std=0.1)).cuda()\n",
    "        fake,_ = netG.forward(noise)\n",
    "        label.fill_(real_label) \n",
    "        output,_ = netD.forward(fake)\n",
    "        errG = criterion(output, label.float())\n",
    "        errG.backward()\n",
    "        optimizerG.step()\n",
    "        D_G_z2 = output.mean().item()\n",
    "        \n",
    "        \n",
    "\n",
    "    print('[%d/%d][%d/%d] Loss_D: %.4f Loss_G: %.4f D(x): %.4f D(G(z)): %.4f / %.4f' \n",
    "          % (epoch, opt_trn.epochs, i, len(dataloader),\n",
    "             errD.item(), errG.item(), D_x, D_G_z1, D_G_z2), end='')\n",
    "    print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Anomaly Detection"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Define basic settings for inverse mapping"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "class ArgsTest:\n",
    "    workers = 1\n",
    "    batch_size = 1\n",
    "    \n",
    "opt_test=ArgsTest()    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "generator = netG # changing reference variable \n",
    "discriminator = netD # changing reference variable "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Define Test Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define settings for loading data in evaluation mood\n",
    "class TestDataSettings:\n",
    "    \n",
    "    def __init__(self):\n",
    "        self.BASE = 'D:\\\\ResearchDataGtx1060\\\\AnomalyDetectionData\\\\NabDataset\\\\'\n",
    "        self.label_file = 'labels\\\\combined_windows.json'\n",
    "        self.data_file = data_file\n",
    "        self.key = key\n",
    "        self.train = False\n",
    "\n",
    "        \n",
    "test_data_settings = TestDataSettings()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(torch.Size([121, 60, 1]), torch.Size([121]), 121)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# define dataset object and data loader object in evaluation mood for NAB dataset\n",
    "\n",
    "test_dataset = NabDataset(test_data_settings)\n",
    "test_dataloader = torch.utils.data.DataLoader(test_dataset, batch_size=opt_test.batch_size, \n",
    "                                         shuffle=False, num_workers=int(opt_test.workers))\n",
    "\n",
    "test_dataset.x.shape, test_dataset.y.shape, test_dataset.data_len # check the dataset shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Define a function to calculate anomaly score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Lambda = 0.1 according to paper\n",
    "# x is new data, G_z is closely regenerated data\n",
    "\n",
    "def Anomaly_score(x, G_z, Lambda=0.1):\n",
    "    residual_loss = torch.sum(torch.abs(x-G_z)) # Residual Loss\n",
    "    \n",
    "    # x_feature is a rich intermediate feature representation for real data x\n",
    "    output, x_feature = discriminator(x.to(device)) \n",
    "    # G_z_feature is a rich intermediate feature representation for fake data G(z)\n",
    "    output, G_z_feature = discriminator(G_z.to(device)) \n",
    "    \n",
    "    discrimination_loss = torch.sum(torch.abs(x_feature-G_z_feature)) # Discrimination loss\n",
    "    \n",
    "    total_loss = (1-Lambda)*residual_loss.to(device) + Lambda*discrimination_loss\n",
    "    return total_loss"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Inverse mapping to latent space and reconstruction of data for estimating anomaly score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 tensor([0.])\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-18-bb890af63b2b>:6: UserWarning: nn.init.normal is now deprecated in favor of nn.init.normal_.\n",
      "  z = Variable(init.normal(torch.zeros(opt_test.batch_size,\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "~~~~~~~~loss=62.367950439453125,  y=tensor([0.]) ~~~~~~~~~~\n",
      "1 tensor([0.])\n",
      "~~~~~~~~loss=138.82664489746094,  y=tensor([0.]) ~~~~~~~~~~\n",
      "2 tensor([0.])\n",
      "~~~~~~~~loss=150.25955200195312,  y=tensor([0.]) ~~~~~~~~~~\n",
      "3 tensor([0.])\n",
      "~~~~~~~~loss=63.134788513183594,  y=tensor([0.]) ~~~~~~~~~~\n",
      "4 tensor([0.])\n",
      "~~~~~~~~loss=38.67523193359375,  y=tensor([0.]) ~~~~~~~~~~\n",
      "5 tensor([0.])\n",
      "~~~~~~~~loss=90.27053833007812,  y=tensor([0.]) ~~~~~~~~~~\n",
      "6 tensor([0.])\n",
      "~~~~~~~~loss=86.46223449707031,  y=tensor([0.]) ~~~~~~~~~~\n",
      "7 tensor([0.])\n",
      "~~~~~~~~loss=159.2848358154297,  y=tensor([0.]) ~~~~~~~~~~\n",
      "8 tensor([0.])\n",
      "~~~~~~~~loss=71.58711242675781,  y=tensor([0.]) ~~~~~~~~~~\n",
      "9 tensor([0.])\n",
      "~~~~~~~~loss=102.59327697753906,  y=tensor([0.]) ~~~~~~~~~~\n",
      "10 tensor([0.])\n",
      "~~~~~~~~loss=108.88626098632812,  y=tensor([0.]) ~~~~~~~~~~\n",
      "11 tensor([0.])\n",
      "~~~~~~~~loss=648.0296020507812,  y=tensor([0.]) ~~~~~~~~~~\n",
      "12 tensor([0.])\n",
      "~~~~~~~~loss=159.79354858398438,  y=tensor([0.]) ~~~~~~~~~~\n",
      "13 tensor([0.])\n",
      "~~~~~~~~loss=85.76463317871094,  y=tensor([0.]) ~~~~~~~~~~\n",
      "14 tensor([0.])\n",
      "~~~~~~~~loss=52.61214065551758,  y=tensor([0.]) ~~~~~~~~~~\n",
      "15 tensor([0.])\n",
      "~~~~~~~~loss=149.47415161132812,  y=tensor([0.]) ~~~~~~~~~~\n",
      "16 tensor([0.])\n",
      "~~~~~~~~loss=69.43896484375,  y=tensor([0.]) ~~~~~~~~~~\n",
      "17 tensor([0.])\n",
      "~~~~~~~~loss=54.840789794921875,  y=tensor([0.]) ~~~~~~~~~~\n",
      "18 tensor([0.])\n",
      "~~~~~~~~loss=152.26856994628906,  y=tensor([0.]) ~~~~~~~~~~\n",
      "19 tensor([0.])\n",
      "~~~~~~~~loss=36.78879165649414,  y=tensor([0.]) ~~~~~~~~~~\n",
      "20 tensor([0.])\n",
      "~~~~~~~~loss=149.66575622558594,  y=tensor([0.]) ~~~~~~~~~~\n",
      "21 tensor([0.])\n",
      "~~~~~~~~loss=142.82687377929688,  y=tensor([0.]) ~~~~~~~~~~\n",
      "22 tensor([0.])\n",
      "~~~~~~~~loss=17.921302795410156,  y=tensor([0.]) ~~~~~~~~~~\n",
      "23 tensor([0.])\n",
      "~~~~~~~~loss=47.16326904296875,  y=tensor([0.]) ~~~~~~~~~~\n",
      "24 tensor([0.])\n",
      "~~~~~~~~loss=88.36796569824219,  y=tensor([0.]) ~~~~~~~~~~\n",
      "25 tensor([0.])\n",
      "~~~~~~~~loss=59.892822265625,  y=tensor([0.]) ~~~~~~~~~~\n",
      "26 tensor([0.])\n",
      "~~~~~~~~loss=114.33897399902344,  y=tensor([0.]) ~~~~~~~~~~\n",
      "27 tensor([0.])\n",
      "~~~~~~~~loss=36.24451446533203,  y=tensor([0.]) ~~~~~~~~~~\n",
      "28 tensor([0.])\n",
      "~~~~~~~~loss=55.80976486206055,  y=tensor([0.]) ~~~~~~~~~~\n",
      "29 tensor([0.])\n",
      "~~~~~~~~loss=104.2334213256836,  y=tensor([0.]) ~~~~~~~~~~\n",
      "30 tensor([0.])\n",
      "~~~~~~~~loss=712.2062377929688,  y=tensor([0.]) ~~~~~~~~~~\n",
      "31 tensor([0.])\n",
      "~~~~~~~~loss=695.2733764648438,  y=tensor([0.]) ~~~~~~~~~~\n",
      "32 tensor([0.])\n",
      "~~~~~~~~loss=696.0294799804688,  y=tensor([0.]) ~~~~~~~~~~\n",
      "33 tensor([0.])\n",
      "~~~~~~~~loss=113.18748474121094,  y=tensor([0.]) ~~~~~~~~~~\n",
      "34 tensor([0.])\n",
      "~~~~~~~~loss=658.4676513671875,  y=tensor([0.]) ~~~~~~~~~~\n",
      "35 tensor([0.])\n",
      "~~~~~~~~loss=94.52996826171875,  y=tensor([0.]) ~~~~~~~~~~\n",
      "36 tensor([0.])\n",
      "~~~~~~~~loss=96.159912109375,  y=tensor([0.]) ~~~~~~~~~~\n",
      "37 tensor([0.])\n",
      "~~~~~~~~loss=102.61505126953125,  y=tensor([0.]) ~~~~~~~~~~\n",
      "38 tensor([0.])\n",
      "~~~~~~~~loss=671.7782592773438,  y=tensor([0.]) ~~~~~~~~~~\n",
      "39 tensor([0.])\n",
      "~~~~~~~~loss=89.71293640136719,  y=tensor([0.]) ~~~~~~~~~~\n",
      "40 tensor([0.])\n",
      "~~~~~~~~loss=126.84182739257812,  y=tensor([0.]) ~~~~~~~~~~\n",
      "41 tensor([0.])\n",
      "~~~~~~~~loss=688.8107299804688,  y=tensor([0.]) ~~~~~~~~~~\n",
      "42 tensor([0.])\n",
      "~~~~~~~~loss=687.2042236328125,  y=tensor([0.]) ~~~~~~~~~~\n",
      "43 tensor([0.])\n",
      "~~~~~~~~loss=691.825927734375,  y=tensor([0.]) ~~~~~~~~~~\n",
      "44 tensor([0.])\n",
      "~~~~~~~~loss=677.1293334960938,  y=tensor([0.]) ~~~~~~~~~~\n",
      "45 tensor([0.])\n",
      "~~~~~~~~loss=114.5341796875,  y=tensor([0.]) ~~~~~~~~~~\n",
      "46 tensor([0.])\n",
      "~~~~~~~~loss=670.5001831054688,  y=tensor([0.]) ~~~~~~~~~~\n",
      "47 tensor([0.])\n",
      "~~~~~~~~loss=88.31558227539062,  y=tensor([0.]) ~~~~~~~~~~\n",
      "48 tensor([0.])\n",
      "~~~~~~~~loss=99.50865173339844,  y=tensor([0.]) ~~~~~~~~~~\n",
      "49 tensor([0.])\n",
      "~~~~~~~~loss=141.00955200195312,  y=tensor([0.]) ~~~~~~~~~~\n",
      "50 tensor([0.])\n",
      "~~~~~~~~loss=677.3421020507812,  y=tensor([0.]) ~~~~~~~~~~\n",
      "51 tensor([0.])\n",
      "~~~~~~~~loss=146.40463256835938,  y=tensor([0.]) ~~~~~~~~~~\n",
      "52 tensor([0.])\n",
      "~~~~~~~~loss=715.7846069335938,  y=tensor([0.]) ~~~~~~~~~~\n",
      "53 tensor([0.])\n",
      "~~~~~~~~loss=705.191162109375,  y=tensor([0.]) ~~~~~~~~~~\n",
      "54 tensor([0.])\n",
      "~~~~~~~~loss=683.5562133789062,  y=tensor([0.]) ~~~~~~~~~~\n",
      "55 tensor([0.])\n",
      "~~~~~~~~loss=690.9383544921875,  y=tensor([0.]) ~~~~~~~~~~\n",
      "56 tensor([0.])\n",
      "~~~~~~~~loss=679.0887451171875,  y=tensor([0.]) ~~~~~~~~~~\n",
      "57 tensor([0.])\n",
      "~~~~~~~~loss=116.42013549804688,  y=tensor([0.]) ~~~~~~~~~~\n",
      "58 tensor([0.])\n",
      "~~~~~~~~loss=691.9901733398438,  y=tensor([0.]) ~~~~~~~~~~\n",
      "59 tensor([1.])\n",
      "~~~~~~~~loss=711.6446533203125,  y=tensor([1.]) ~~~~~~~~~~\n",
      "60 tensor([1.])\n",
      "~~~~~~~~loss=696.644775390625,  y=tensor([1.]) ~~~~~~~~~~\n",
      "61 tensor([1.])\n",
      "~~~~~~~~loss=744.7349243164062,  y=tensor([1.]) ~~~~~~~~~~\n",
      "62 tensor([1.])\n",
      "~~~~~~~~loss=771.4285888671875,  y=tensor([1.]) ~~~~~~~~~~\n",
      "63 tensor([1.])\n",
      "~~~~~~~~loss=713.0824584960938,  y=tensor([1.]) ~~~~~~~~~~\n",
      "64 tensor([1.])\n",
      "~~~~~~~~loss=694.4222412109375,  y=tensor([1.]) ~~~~~~~~~~\n",
      "65 tensor([1.])\n",
      "~~~~~~~~loss=708.8585205078125,  y=tensor([1.]) ~~~~~~~~~~\n",
      "66 tensor([0.])\n",
      "~~~~~~~~loss=696.4937133789062,  y=tensor([0.]) ~~~~~~~~~~\n",
      "67 tensor([0.])\n",
      "~~~~~~~~loss=120.4222412109375,  y=tensor([0.]) ~~~~~~~~~~\n",
      "68 tensor([0.])\n",
      "~~~~~~~~loss=682.6812133789062,  y=tensor([0.]) ~~~~~~~~~~\n",
      "69 tensor([0.])\n",
      "~~~~~~~~loss=135.65298461914062,  y=tensor([0.]) ~~~~~~~~~~\n",
      "70 tensor([0.])\n",
      "~~~~~~~~loss=724.67236328125,  y=tensor([0.]) ~~~~~~~~~~\n",
      "71 tensor([0.])\n",
      "~~~~~~~~loss=688.5003051757812,  y=tensor([0.]) ~~~~~~~~~~\n",
      "72 tensor([0.])\n",
      "~~~~~~~~loss=666.705810546875,  y=tensor([0.]) ~~~~~~~~~~\n",
      "73 tensor([0.])\n",
      "~~~~~~~~loss=125.47927856445312,  y=tensor([0.]) ~~~~~~~~~~\n",
      "74 tensor([0.])\n",
      "~~~~~~~~loss=667.2129516601562,  y=tensor([0.]) ~~~~~~~~~~\n",
      "75 tensor([0.])\n",
      "~~~~~~~~loss=92.22420501708984,  y=tensor([0.]) ~~~~~~~~~~\n",
      "76 tensor([0.])\n",
      "~~~~~~~~loss=88.28347778320312,  y=tensor([0.]) ~~~~~~~~~~\n",
      "77 tensor([0.])\n",
      "~~~~~~~~loss=73.57276916503906,  y=tensor([0.]) ~~~~~~~~~~\n",
      "78 tensor([0.])\n",
      "~~~~~~~~loss=631.4630737304688,  y=tensor([0.]) ~~~~~~~~~~\n",
      "79 tensor([0.])\n",
      "~~~~~~~~loss=96.48371124267578,  y=tensor([0.]) ~~~~~~~~~~\n",
      "80 tensor([0.])\n",
      "~~~~~~~~loss=111.462158203125,  y=tensor([0.]) ~~~~~~~~~~\n",
      "81 tensor([0.])\n",
      "~~~~~~~~loss=79.0983657836914,  y=tensor([0.]) ~~~~~~~~~~\n",
      "82 tensor([0.])\n",
      "~~~~~~~~loss=102.3183364868164,  y=tensor([0.]) ~~~~~~~~~~\n",
      "83 tensor([0.])\n",
      "~~~~~~~~loss=86.724609375,  y=tensor([0.]) ~~~~~~~~~~\n",
      "84 tensor([0.])\n",
      "~~~~~~~~loss=33.38389205932617,  y=tensor([0.]) ~~~~~~~~~~\n",
      "85 tensor([0.])\n",
      "~~~~~~~~loss=61.75398635864258,  y=tensor([0.]) ~~~~~~~~~~\n",
      "86 tensor([0.])\n",
      "~~~~~~~~loss=82.82698822021484,  y=tensor([0.]) ~~~~~~~~~~\n",
      "87 tensor([0.])\n",
      "~~~~~~~~loss=80.24472045898438,  y=tensor([0.]) ~~~~~~~~~~\n",
      "88 tensor([0.])\n",
      "~~~~~~~~loss=645.3245239257812,  y=tensor([0.]) ~~~~~~~~~~\n",
      "89 tensor([0.])\n",
      "~~~~~~~~loss=49.434974670410156,  y=tensor([0.]) ~~~~~~~~~~\n",
      "90 tensor([0.])\n",
      "~~~~~~~~loss=148.80616760253906,  y=tensor([0.]) ~~~~~~~~~~\n",
      "91 tensor([0.])\n",
      "~~~~~~~~loss=146.02944946289062,  y=tensor([0.]) ~~~~~~~~~~\n",
      "92 tensor([0.])\n",
      "~~~~~~~~loss=166.4810791015625,  y=tensor([0.]) ~~~~~~~~~~\n",
      "93 tensor([0.])\n",
      "~~~~~~~~loss=57.47191619873047,  y=tensor([0.]) ~~~~~~~~~~\n",
      "94 tensor([0.])\n",
      "~~~~~~~~loss=31.61819839477539,  y=tensor([0.]) ~~~~~~~~~~\n",
      "95 tensor([0.])\n",
      "~~~~~~~~loss=169.37582397460938,  y=tensor([0.]) ~~~~~~~~~~\n",
      "96 tensor([0.])\n",
      "~~~~~~~~loss=50.99140930175781,  y=tensor([0.]) ~~~~~~~~~~\n",
      "97 tensor([0.])\n",
      "~~~~~~~~loss=51.639198303222656,  y=tensor([0.]) ~~~~~~~~~~\n",
      "98 tensor([0.])\n",
      "~~~~~~~~loss=155.75277709960938,  y=tensor([0.]) ~~~~~~~~~~\n",
      "99 tensor([1.])\n",
      "~~~~~~~~loss=39.53273010253906,  y=tensor([1.]) ~~~~~~~~~~\n",
      "100 tensor([1.])\n",
      "~~~~~~~~loss=155.7327117919922,  y=tensor([1.]) ~~~~~~~~~~\n",
      "101 tensor([1.])\n",
      "~~~~~~~~loss=158.93218994140625,  y=tensor([1.]) ~~~~~~~~~~\n",
      "102 tensor([1.])\n",
      "~~~~~~~~loss=164.66543579101562,  y=tensor([1.]) ~~~~~~~~~~\n",
      "103 tensor([1.])\n",
      "~~~~~~~~loss=210.72886657714844,  y=tensor([1.]) ~~~~~~~~~~\n",
      "104 tensor([1.])\n",
      "~~~~~~~~loss=149.5918731689453,  y=tensor([1.]) ~~~~~~~~~~\n",
      "105 tensor([1.])\n",
      "~~~~~~~~loss=53.95457077026367,  y=tensor([1.]) ~~~~~~~~~~\n",
      "106 tensor([1.])\n",
      "~~~~~~~~loss=172.0368194580078,  y=tensor([1.]) ~~~~~~~~~~\n",
      "107 tensor([0.])\n",
      "~~~~~~~~loss=39.402374267578125,  y=tensor([0.]) ~~~~~~~~~~\n",
      "108 tensor([0.])\n",
      "~~~~~~~~loss=30.838376998901367,  y=tensor([0.]) ~~~~~~~~~~\n",
      "109 tensor([0.])\n",
      "~~~~~~~~loss=138.29322814941406,  y=tensor([0.]) ~~~~~~~~~~\n",
      "110 tensor([0.])\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "~~~~~~~~loss=144.58026123046875,  y=tensor([0.]) ~~~~~~~~~~\n",
      "111 tensor([0.])\n",
      "~~~~~~~~loss=177.72926330566406,  y=tensor([0.]) ~~~~~~~~~~\n",
      "112 tensor([0.])\n",
      "~~~~~~~~loss=155.30294799804688,  y=tensor([0.]) ~~~~~~~~~~\n",
      "113 tensor([0.])\n",
      "~~~~~~~~loss=45.42596435546875,  y=tensor([0.]) ~~~~~~~~~~\n",
      "114 tensor([0.])\n",
      "~~~~~~~~loss=169.12242126464844,  y=tensor([0.]) ~~~~~~~~~~\n",
      "115 tensor([0.])\n",
      "~~~~~~~~loss=38.48234558105469,  y=tensor([0.]) ~~~~~~~~~~\n",
      "116 tensor([0.])\n",
      "~~~~~~~~loss=46.81413269042969,  y=tensor([0.]) ~~~~~~~~~~\n",
      "117 tensor([0.])\n",
      "~~~~~~~~loss=217.19903564453125,  y=tensor([0.]) ~~~~~~~~~~\n",
      "118 tensor([0.])\n",
      "~~~~~~~~loss=163.85986328125,  y=tensor([0.]) ~~~~~~~~~~\n",
      "119 tensor([0.])\n",
      "~~~~~~~~loss=77.74629211425781,  y=tensor([0.]) ~~~~~~~~~~\n",
      "120 tensor([0.])\n",
      "~~~~~~~~loss=178.9587860107422,  y=tensor([0.]) ~~~~~~~~~~\n"
     ]
    }
   ],
   "source": [
    "loss_list = []\n",
    "#y_list = []\n",
    "for i, (x,y) in enumerate(test_dataloader):\n",
    "    print(i, y)\n",
    "    \n",
    "    z = Variable(init.normal(torch.zeros(opt_test.batch_size,\n",
    "                                     test_dataset.window_length, \n",
    "                                     test_dataset.n_feature),mean=0,std=0.1),requires_grad=True)\n",
    "    #z = x\n",
    "    z_optimizer = torch.optim.Adam([z],lr=1e-2)\n",
    "    \n",
    "    loss = None\n",
    "    for j in range(50): # set your interation range\n",
    "        gen_fake,_ = generator(z.cuda())\n",
    "        loss = Anomaly_score(Variable(x).cuda(), gen_fake)\n",
    "        loss.backward()\n",
    "        z_optimizer.step()\n",
    "\n",
    "    loss_list.append(loss) # Store the loss from the final iteration\n",
    "    #y_list.append(y) # Store the corresponding anomaly label\n",
    "    print('~~~~~~~~loss={},  y={} ~~~~~~~~~~'.format(loss, y))\n",
    "    #break"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualise Anomaly Detection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD/CAYAAADoiI2GAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABsXklEQVR4nO19eZgkVZX9eRGZVVn73tVL9d4NTdPdNNAgyCKCgCKyCC6M4yCLzE8HZXBGB5cZx1FHBWWcGRTHjUURHBRGRTZFsdmhgabpjd6X6qW6qrv2NTPi/f6IfBEvXrxYcqvMrH7n+/iozsrMeJEVeeO8c8+9l1BKoaCgoKBQftCKvQAFBQUFheygAriCgoJCmUIFcAUFBYUyhQrgCgoKCmUKFcAVFBQUyhSxyTxYa2srnTdv3mQeUkFBQaHs8eqrr/ZQStvExyc1gM+bNw9r1qyZzEMqKCgolD0IIbtljysJRUFBQaFMoQK4goKCQplCBXAFBQWFMoUK4AoKCgplChXAFRQUFMoUKoArKCgolClUAFdQUFAoU6gArqAgwaHBMXzj0U0YTxnFXoqCgi9UAFdQkOCnz+7C/6zegVd39xZ7KQoKvlABXEFBgGFSPPx6JwBga9dQkVejoOAPFcAVFAQ8t60HXQPjAIC3ugaLvBoFBX9Mai8UBYVywK9f60R9IoZ5rTXYqgK4QglDMXAFBQ6DY0k8seEg3nfCTCyb1YAtXUNQc2MVShUqgCsocHjszYMYS5q44uQOHDOtFv2jSRwaHC/2shQUpFASioICh4de78SC1hqcOLsRY0nLQrilaxDt9Ykir0xBwQvFwBUU0hhPGXh1dy/OX9oOQgiOaa8DALx1UOngCqUJFcAVFNJ46+AgkgbFCbMbAQCttZVoqalQVkKFkkV5SCiP3QIcfLPYq1CY4pg2MIYHKoax8sVGYI0OAPiZ3g/zLQB3NRR1bQpTANOXA+/5Zl7fMpSBE0J+Sgg5RAhZzz12GyFkMyFkHSHkYUJIY15XpaBQBAyPpxDTCCpjzteiqiKG0QkDFMqJolB6iMLA7wZwB4B7ucf+AODzlNIUIeRbAD4P4J/yv7w08nzXUlCQ4e++uxrT2xK4+5pT7ceee3E3vvR/6/Hc5ediVmNVTu//3T9uwcUrZmDRtLpcl6qgACACA6eUrgZwRHjsSUppKv3PFwF0FGBtCgqThpGJFLZ0DWLFLLdUwhKZW3JMZA6Pp/DdP27Fj5/ZmdP7KCjwyEcS81oAj/n9khByAyFkDSFkTXd3dx4Op3C0Y/fhYWzY35/X99ywfwAmBVZ0NLoeP6a9FoBlJcwFKdOSYJ7Z2qMKgxTyhpwCOCHkiwBSAO7zew6l9IeU0lWU0lVtbW25HE4hIg4NjOGJDQeLvYyC4ZuPbcbl338ez2/vyel9dnQPYWjc2kiu67RuCCs63Ay8sboC0+oqsTlHBm6kA/i+vlHs6BnO6b0UFBiyDuCEkKsBXAzgI1RRipLCL1/Zi7/92at4bltuAa5UMTJhYCJl4uP3rMEbe/uyeo+D/WN4z38+g2vuehmmSbGusw/T6xOYJinYOWNRKx5ffxCHBsfsx57a1IXfvbE/8vFSpmn//MwWtRNVyA+yCuCEkHfDSlpeQikdye+SFHLFeMoKFl/53QYkDTPk2eUHw6RY0FqDltpKXH3Xy9h9OHNGe+fT2zCeMvHKrl7c+8IurOvs97BvhpvOW4ykYeKOP20DYMkpn7jvNXz6gdfxp81dkY6XMhyO88zWqXljVZh8RLER3g/gBQDHEkI6CSHXwXKl1AH4AyFkLSHkBwVep0IGYHrrlq4h/PzF3UVeTf6RNEy01lXi59e9DWNJww6sUXGgfxT3v7wXHz5lNs45tg3fevwt7OwZtgt4RMxrrcEHT5mNX7y0B9u7h/CZ/12LusoYjptej5vuX4sd3eGFPkxCqUvE8MKOw5hITb0bq8LkI4oL5SpK6QxKaZxS2kEp/QmldBGldDaldGX6v/83GYtViAbDNFEV13HW4lbc/octeGXXEaze0j1lpssYJkVMI5jTUo0PrpqN/1u7D10DY+EvTOP7f94Ok1L83TsX4d8vXw5dIwCA5bP8i3VuOm8xdI3ggz94Aev3DeDrly/Hj65ehXhMww0/exXD4ynf1wLOTfXsxW0YmTDw2p6p8bdQKC5UKf0URMqkiOkEX37fUoxOGPjAD17A3/z0ZXzgB89jYCw5KWsYHEti75HCqGtJkyKmW5fu9WcugGFS3PXcLvv3QSmZ/X2j+OUre/GBVbMxu7kaMxur8JVLjseMhgRWzmn0fV17fQLXnDEfh4cn8P6TZuHdy6ZjVmMVvnXFCmw7NBSab0ilpayzFrdC1wie2ap0cIXcUR6l9AoZgTHURdPq8Mu/PR2Hh8bx8s4j+PGzOzE2YaA+Ec/r8XYfHsbTb3VjXmsNZjQk8PDr+/DzF3djPGnikU+faXup8wXDNBFLs+Y5LdV4z7IZuO+l3bjx3EV4alMXvvK7jfjkOQtx/VkLPK99ZN1+TBgmPnnOQvuxK07uwBUnh5cy3HjuIrTWVuBDp8y2H5vXUg0AmAjJNTAG3lgdx0lzGvHM1h589sLwc1VQCIIK4FMQSYNC1yyGevLcJgBA78iE9Tsz/4ahHz+zEz/jtHaNAO9ZNgMv7jiMf/jfN/DQJ9+OuJ6/zV7KoHYAB4Abzl6A3795AFfe+Tw2HxxEfSKGrz+6CXNbanD+0nbXa0cnrECbTVVlbWXMc1NgOwEj5HNlv9c1DeccOw23PfEWtncPYWFbbcbrUFBgUBLKFATPUBlYQDeM/Afw8ZSB1tpKPHDDabj1ihV46h/Owfc+chK+dtkyvLmvHz94entej8ckIoYTZjfitAXN2HxwEJ84ZyGeu+VcLJ/VgJseeB2bDgwIrzVBCKAJn0+2YJ9zkvtcx5IGPvg/L+DNTqfYiLmBYjrBB1fNRkVMw0+eVVWZCrlBBfApiJRJ7cQcQzwd8Hg/ct6OZ1Ak4hpOW9CCD54yG/NbawAA71k+A5ecMBP/9aetnkAaFWt2HfEk/CyJyH3p/vdVJ+GRT52Jf3r3EtQl4vjR36xCXSKGT93/unutJkVcy99lz24kBve5dg9aktW6fX2uNQNWwG+rq8QVJ3XgV692omdITftRyB4qgE9BGCa1AzYDC+ipLCWUVIDGmzKpr0TylUuOh2FS/H7dgayOe+vjb+H2J7e4Hksa3h1GW10llnEukvb6BC47cZYnkZoyTM/NLRfoEgbOgjXv/U7ZEor1/OvPmo+JlIl7X5h6Nk+FyYMK4FMQMgbOAl4qCwllf98olv7LE1i/T95/JGX6B8WmmgpUxDRPQdF3nnwLf958KPTY4ynD45k2BAnFD3FN89ywRPklV7CdAK+Bs10Of2z2ubPnL2yrxbuOa8fPXtiF0Qkjb+tROLqgAvgUhGF4JQb272wklIMDY5gwTHT2jkp/LyYVRcQ0zcVQAeDnL+7GkxvDqxiTBkVSWDOfpA2CrhEYJnXZCplDJ1+I6d6dTcpm4Cb3mOl6PmAlX3tHknj49X2u97z5l2vxi5f25G2NClMXKoBPQcgYsS4JNFHB2KWf0yKM1cZ04tKIASvoi4/5HVvcNciStDLEJeccNfhHhbOz4YJ1er38cXkNnOGUeU2oiuvY2eOu5Hxma48q9FGIBBXApyBkAZUFDj4IHxocw7/+NrxfCvu9H3tPSZKK4rFF+2JKEpilxzZNrwxiRJNBdIm8YZimJz+QC5ydjYyB+2vgAEAIQUwnEpnHDLUlKigAKoBPSRhSDdz6U/PB+vlth3H387uwK6S9qSwpxyMlSSqKxxbtiylJYJa/N/UkUFMRZRAWqJMCO85nElOWW0hJbnjs92KyN6YRz+eaMuiUbEKmkH+oQp4piJThtco5djdeTvAm26Tvl6OEomvEpWNTSpE0aCSWKXONpEzTLqAJgi7ZdQQ5ZrKBphEQ4rYR2gxcktj03Fh1zbOzUQxcISoUA5fgtic246YHXg9/YolCxsBlNsIwZs0g03TdvzcDJZS4TlwBif0YhWWmTC8bjcrAWZAXLX75ZOCA1+1if17cumUauPVaOQPP1u6ZbyQN01WQpFBaUAFcgk0HBrFxf3aFJ6UAi6GKgSKt1XLBImkzxeBAapjBGniYrU8XghQLyJEYuKCVW64SBN4wGGS6v8xDnit0jUjZdtLwBnVx56ALGjilNH3OpSGhPL7+IC793rOq4KhEoQK4BCmzdBhQNpD5wB05gddlo0koSSOYqSdDbIRxQSZgATVKX5akYUqDY5QkplPm7j52Pn3gbC0pSbAWpRt+TQxxTfOsj39+sdE/moRJgZFx5VUvRSgNXIKUYZZ1Eknmy5b5laNKKGE2QllpOw+RgTsBLoKEYlDoWrgUIYNU9zfzayNka0lJNXB+3X4aOJEG+mwKrgoB52Zbvt+HqQzFwCUodxeA3IXidUvYzDrky8k+C78vcVIi2biOrWtSFh0lSIk+cLbmKDq2rHjJME3E8y6hhJ8fW7cnuSwUOWUiL93/8p6CF/xksh6FyYcK4BKkTLNkGFA2kLk0pH7liBKKzcB9bYRhlZh+DDWaD1wmMURxksQkidt82wittRDXZyOTQex2sqI/XydSeSlKxezDr+3Dw693Zr/wCIi6S1MoDlQAl0DmfCgnyMrFdUnXvKjb9WRIwLV05eBCHllRS5QbB6WClmzIpQjpcXVv4jbfvVDYWvjdCWPULv+5j/QT04iQZM305lbYwBpmIVUoLlQAlyBZQjaubCBtJyvpmsdYXpgWbYRUYoY5O8Rqw5QR7bh8UyjWz8QvEEqPK2PgIXp9Nojrmqfa0/p/+I0npmtCoM9QXiqwNs3WoTTw0oQK4BIY5S6hSCQNv6IWAKEsLowxhzk7YoJGHOZqsY8rY+0+djzpcVni1lWJWSAboUTvTkrWL2PgsgRvlMCcNKK1I8gFKcnNSKF0oAK4BCmDYsIwA4fjljJSEqeFXAOPtj22t9G+NsLgQh4rSGVulQsKbFGCsF/xUt5thFq4js0Sy4RIGLjMhRIhYBoR2xHkglJzxSi4ERrACSE/JYQcIoSs5x5rJoT8gRCyNf3/psIuc3KRLHPWIevW58dGgfCKyLCAG9aiVbTKRXU2uHRlocd2pH7gEg087GaTDaKcn0zWAixpS+rNjxAwZX1i8o2UcqGUNKJcyXcDeLfw2C0AnqKULgbwVPrfUwZhpeOlDlmiTsZGoyaowroRJk3qcVfwiPkUq0S9cQAO+3cGI2TCwL1MOJ/QBSug7PxShty+6PHIZ5A0nMwkptLASxOhAZxSuhrAEeHhSwHck/75HgCX5XdZxUUqYoApVcgYcVwyPd1JUEWzEfqxQsMMnjPpLVaJyMANLwN3CnnCuYesfUAhXChxLbwYx5eBi0nMDDRww5i8JGYhhmEr5I5s95LtlNIDAJD+/zS/JxJCbiCErCGErOnu7s7ycJOLTLaxpQiZBs5iR0oWFEP7gfvvSCiloaxW14i0L0g2GjhbcxDj548rHifMs54NrPMLH+ggS7x65JcMXChJM1pHx1yQiSavMPkoeBKTUvpDSukqSumqtra2Qh8uL7CDRZkycJnTghCSTrZ5t/rhfmx/xsxeGzQkIa5pPgw1mo3Qeq5bWogyWd6ZyOP2Y0dxsGQCvx2G+7hy94t4c8ukF4ph0oJLKEF/e4XiI9sruYsQMgMA0v8Pn05bRnB0v/K7aE2TwqTyQhevHztasJAVpojvEdRfRNezq8SUd/iLXsgjt07m30Yo2iSlEooP87da0fJFQIyBR7ERmpOQxIwu6ShMPrIN4L8FcHX656sB/CY/yykNOMUU5XfRGtQ/yRfTNGlb17DzDGpmxT6rYAbuV8iTgYQirCETF4qL4RaglN7TKkDSosDwSfR62HsGCXTDpAUnGcpGWNqIYiO8H8ALAI4lhHQSQq4D8E0A5xNCtgI4P/3vKQE2LQYIL3ApRTgBzvun1QXLWmQ/tukNSPbvIjSX0j03jmhBwZXEFPISmbhQxPYB+ZzIA0jayUokoqRPBWjMR0KJNq1oMjRwJaGUMkLbyVJKr/L51Xl5XktJgL9Oy3HbGFRqHtfdw4WTEdmVbMKM53gBQTGeZcOmIGYaaaCDPRPTLaHkn4H7TORxMXC5dBPTNXeg59g7pdRT+MMjZZowKUKflwvK3VI71aEqMQW4WF+q/C5aI4AR65q7a17UboRBmnWUykiv19l6jUktzd4Psr4tbM3RBjpIrJMmLUA72Qg2Qh/pxpOXkCSZZWC5DqCwO0VZb3OF0oEK4AJcPTvK8KJla5YFOP9kW5gbJEADjyBpsH7gdkMqidwgP66/NS+TgQ58ZSSlwQnXbBDTRRuhzIUi95+LziC/n0XIEryFgNLASxsqgAtISfy85QS777Qv2wtOtsnf09+bHKW0XZxNKZNTZJC3oGU3qMxnYmYyji0TiC1hZbkFvy6IsbTF0rm5ed06MshcPYWAKqUvbagALkDmlignBDFU7/DdaNvjoMk9hi2hBFdiuo8XbZfj7mCYeRJTbOCVyTi2TCBOHJLlFvw08Lig07teE0AgZJp7IaAKeUobKoALkE1tLycYdoCTOx4Mw/vFD2NXQc6IZBQJRaiIjBykJLuhTGyE4hi5TMaxZQKx26JUQvHVwN06feSbmySXUQikMvClK0w+VAAX4E5ilt9FGyQTxISiEac4JpqdT/a8INsif1yAa0gVUeeVtVllj0UJwppGoBGvFS7fNkLfnY0gAcmOy24yTrfFqPLS5EgomVSGKkw+VAAXICuJLieEa+DeQJMbAw93odjJREmBVNBnLHse69sSpZTeWpfGBdToVZyZIK5rUkYsThKS/k18dgnWz9HkpUJKKMmIuzSF4kAFcAHusubyu2iDJA2P3S1iz5egcuookoZo54safIK81VGaWbF1icnaQjSzktkIAf7m5+8Dt15j2s8TXyuDW+orHNFQDLy0oQK4AJn3uJzgMHDvnzauia1LI5a0B1RiRtGVbZlA0uUxcpDKMgjz8kYUuScbxD2l9HI5JVMGnqnFshCQ7ZoUSgcqgAuQTQgvJwQV1vgxxcgzMQM08CBdmbFzmY0wUEKRNcCyJZtoly4vb0SRe7KBrmmuoiTZdB5fDVyYGiRroSvDZPnAFQMvbagALkDWf6OcECRp+Gvg2Usodn/ukEpMQN5XPJBlSn3g+WDg+bYRul02SYkM4tczXWx5a0S9uUUM9LkiqlNJoThQAVzAZCWHCgVbI/bZrgfJEv7v6S+1MGdJUFJRnAYUtVhKWnRkUGjEcphEQZyz+EWxPGYDxybp/Zz4Vrx+uyLrtW6XDZDBZ1PQSkyvO0ahdKACuABZgCsnBI0c04VS+rBRaQxReqFEY+DeG0agjdBnik8mQ4l1rl1rUH4gF3iCsCSPYk3kkVs7rdfIHDrRJJSC9kKJ2DlSoThQAVzAZJUoFwpBAw/iujgB3V8a4REU6CNN5BE1cKFC0Q8yp4VhmhlJIHFNcyojC1hKDzi7Efln7B1zB0g+G8kNVoaoieBcEdVqqlAcqAAuYLIq3AqFoHJxsSugUzwSwsADAn3UfuD866PKVC5LHefSyMTHzfdAL5SNkCUiZfMs+QAYJKHISumDrj/ZFJ9CIGq/HIXiQAVwAfyXYaIMt42hGnj696ZpdeYDokgo/hq4w8ADNHDBKhdVJkiaFBV2cHQCYSaVlDFd8wTHTCSYSMfwNM2iHLNm+rt85xC3XSjyDoZ+mLQkpsQ5pFA6UAFcgLvQpfwu2iBbHz88IJmBVMR7lFllofO76Bq4TEsPDlJW0ON7jWQ6kIHvFFgwCcVjBTSRiOkA3JWM8jF3/vJSUNCcrIph1U62tKECuABZ/41yQlQGHlWH5t8TcE8s4n8XWIlpVxtKAnjATTJpWEGPtz/6DQf2PzbxJE8L0cyKf/+USVEZtwJ4mAbutBmQfTZBCd7C52oopZ4EsEJpQQVwAbLxVuUEp71rcCFPJnbJoN4lUQprnGrDzJwWlnNDc1WQ+g1GCDq2yG6j9lGJCnH2ZsqkSMQF3d8wpYleu+Wt4UgtDIFJzEmwu0bdKSkUDyqACyh3H3hQaXtc1zK2q7Hfs5GL4mcSjYF7GaouSAfy41read4KmKmNMMYNVDYiWB6zgbent4kEY+Amtcef+TUY41/LF/wEtTOW9UrPN8q9sdvRABXABUwZH7gkoPoz8HAJhWm6YrCP4uxwWCavEbv9zzIkDSthGdM0V3DMVEIR3S9BlsdsoEuadTEGnjRMGNT/MxIbfSUNan82wRbLwkso7r455UdmjgaoAC6AfdkrY1pBCyQKhcgaOPe8oADAdNDKuHdAMP/vaBN5nEDKNOIgBs7YaFzoKJhpElPU/fOugQuzN1MmRRV3fs5x/TVwvuKRfTZB199kyBuyMXEKpYWcAjgh5GZCyAZCyHpCyP2EkES+FlYssC9bVYVeltvGoH7Z7mSg9bxETAsMAOx3NgMXGHOkfuAeG6HDMoMrMdMuFGHdmdgIdU5CSUW42WQDr5OEk1AMal9HMuYftysxHQmlKsLNbTLqFaJ2RlQoHrK+kgkhswB8GsAqSukyADqAD+drYcUC+7JVxfWy3DYG9cvWuQG6dmCO64EBwLCfJw+4hhnem0QcG2aY1BXgfM/FoIhrmnsoQ4YMPM5LKOxmk3cJxdH4md5dyUlOYUM2ALdfnO12gv8uhZc3lAZe+siVisQAVBFCYgCqAezPfUnFBWMdVXG9vDXwEM8x+9In4no0Bu7DCpMGDe2vLY4NS3IyQaDOa3p94AZXJBMFukQ2ynclZpzzgafEG55hBjbRcnql87sTeb6Bh2tyT4GCq9LASx+xbF9IKd1HCPk2gD0ARgE8SSl9UnweIeQGADcAwJw5c7I93KSBD2zlGMADNXDODcLOrTIeIqHYz2O6rPsz8Zs04zqupFjFj9GL5xLTCEyhmjIzBs67UFiCt1DNrExux5L+vDgGLjtuzFOJ6Xw2YfkB++dJYOBhGngymURnZyfGxsYKspajBYlEAh0dHYjH45Gen3UAJ4Q0AbgUwHwAfQAeJIT8NaX05/zzKKU/BPBDAFi1alXJ38YZIwwLbKWKoHJxvuDEDjSxYAnF0cDlQYUV2wQhJui8KZOitjLmWq/fucR0DRSmK8lXE49+2fLOm6BGX7kgzrlsGBt2grCz9rCbqvV/E/WJeHq9AQyc74VSoOs0amMtAOjs7ERdXR3mzZsHQvL7+R4toJTi8OHD6OzsxPz58yO9Jhcq8i4AOyml3ZTSJICHALw9h/crCSTTW3Rx/Fi5gN2AZDHKtrtxibWqCt01TUYEvyMBvEGFFdsEQdR5rSRfxCSmxiSU7DXwJCe/AIWZicnWxtgwX0ofRdZKuXYn4fISz7oLlcR0ObJCZJqxsTG0tLSo4J0DCCFoaWnJaBeTSwDfA+A0Qkg1sf5q5wHYlMP7lQQsn7GWHoZbhgw8fQOSfZHsghPT0WXDAmnKwyhFH3i4hCJ23OOTmIFBynR84E4iMrNSepn3vVATeQyTZ+BOktbWwGUSiuYvoQQzcOt3lTGtYBY//uYdRaZRwTt3ZPoZZh3AKaUvAfgVgNcAvJl+rx9m+36lAqa7xnWtYFvTQsJvdBcAV/UjL6Gwx2RgX2LHVSGW0ocHVHEiT9I0URmLFqR01guFa2mbaSWmOEy5UDZCXgOvqnBuUJkxcMsmqZEwDdx0rtMCl9InylROPBqQ05VMKf0ypXQJpXQZpfSjlNLxfC2sWLB0V6t4JJkqPwklqNScDxYsqDnJNvm5elwVkkrMMAmFxS3bSZLWtsUhy55jp/uHxHTN1exJZpH0A98LxTBNEFKIZlacC8VwmDFg3aCCNHBNI9CIc3Oxcgqa1TkyxAfueOQLVUrvXCPlUMhTW1tb7CVMOlQlpgDLuubetpcTUoZ/u1V+u24Hmriji0vfTyIJuH8fzsAJSVdT2gzcknnCqkBT6WDGz7W0vOGZlNI7wT8ZYa3ZgK+mlNkuHebv83fRNW64Bq/7h3VqdLcZyDeS3M2oHPNBRwOydqFMVTBJoJw1cP9A4TBwMdD4MnCBUXo0cJ9BBSLE6fC6RlyBWX4uJnSdwKTuuZaZzLR0M/DMEqCZHMNaL3UqXLnP1dHe5euOa4Qbx2Yx67CbGxstJ47Jyyd4S2QmDPwrv9uAjfsH8rqWpTPr8eX3HR/puZRSfO5zn8Njjz0GQgi+9KUv4UMf+hAOHDiAD33oQxgYGEAqlcKdd96Jt7/97bjuuuuwZs0aEEJw7bXX4uabb87r2gsJFcAFsHalcT08816KiKqB8xWn7DEZwgK9X59rEXGupD2ZThSHMnAzzbb5LoqmvC2rH2LpToaUWrJRvlvJAuLnKhbyBLtQ2Ovt3Um6gVdcD94Bst2EOCYvn7BltlhwsVep4aGHHsLatWvxxhtvoKenB6eccgrOPvts/OIXv8CFF16IL37xizAMAyMjI1i7di327duH9evXAwD6+vqKu/gMoQK4APYlj5cxA/frFcJPQHccBu5OgSIMQUIRpRa/PtcidFdJe9qqGRKk3D7w7Ap5RO97Jvp55GPonN6d/nzYKDielfutW2zzq6cDc2ASk8lLBUy2s+NXxjNzukRlyoXCs88+i6uuugq6rqO9vR3veMc78Morr+CUU07Btddei2QyicsuuwwrV67EggULsGPHDnzqU5/Ce9/7XlxwwQVFXXumUBq4ABYgYnp56n5BDDwmY4ohZdtJQUKRJTGjBFS+nwmTQcKCVJIlMYWGVJnOxAScMvd8O1AA/nN1inYsFm1JRGFtbN0um/QOMIRZJ9Oj5cK08lzg/O0tBi6O0ytV+K3z7LPPxurVqzFr1ix89KMfxb333oumpia88cYbOOecc/C9730P119//SSvNjeoAC6AuSriGinbdrJhGrjFFMXkpDwIiKXhMhthFFmCDzRJkwXm4CDFt5MVGWpU8Ba/THuJRz4G/7lyXnMmjTjNrPx3Rny/lrimpXcswZ8Nc+gUup1slNL+UsLZZ5+NX/7ylzAMA93d3Vi9ejVOPfVU7N69G9OmTcPHP/5xXHfddXjttdfQ09MD0zRxxRVX4Ktf/Spee+21Yi8/IygJRQDTWPkBwOWEKC4Uq+AkaiGP11Xh/n00XzZrCWuaFJTC3uWEV2JqoNTxUkdxvbiOy7WyzbSKM/oxnM+V7znOdP+wlrvMCsh6r1vMOtxGWGgG7nEgmRTpDVtJ4/LLL8cLL7yAE044AYQQ3HrrrZg+fTruuece3HbbbYjH46itrcW9996Lffv24ZprroGZPtdvfOMbRV59ZlABXADvQinHQp6gIOVquiQ0qfJjwnYzq5hcK7cqB6NIKMTVK8SqsIyQxNQJgOxnYuouLTqzToZRYfvcTb73t2YH5qApSYDz2fBSS1hgtoiGZt8YCwExT1LqDHxoaAiAZVu97bbbcNttt7l+f/XVV+Pqq6/2vK7cWDcPFcAFMEZpsafyY+BGQIBjwcvd9lQujTCEPS9qaXtM11wVoE6iLriRlp5m4EyDzdRGGOdvWgVi4LbPndO7rfPT3LJKwM7Ikngcu2EUjzxj+YXK1dgMPBZ8k1coHpQGLsCpxNRg0tJnHSKCEnW83c3uhRIyGcfjVsmiEhOwgldS6I0dCykD5ysxXQw1w0Iedh6MtRYCLCHLPi/mZOJL6X01cN29O2El8mFT6aNUs+aCqDd5heJBBXAByfQXQ5xzWC4I6s9tV2Ka1Ga+rGeHr4QSVokZMTHI/NiGLRNoriIbEWyyjVVt6JYiMi2lt86DZmxBzATxNNtmnyvbYUTTwC0roMHd3HSNBF57dsWmXrhke1gRl0LxoSQUAU6TIHeToXJBUJCyKzE5JuzXpIp/P+t5/iPVojBwXbOCFO+JDpIJbDaqE1AQmBQYT/nP+/QD38o20wRoJtDTbJt9ro7PnYZq4PG0lOScM2PvYW0Gwr30ucDLwMvru3A0QDFwAXYSU2jzWS4ISvJ5ilrSDI49JoNoIxQ162SEdrIAH6T4AOefqOOrF5nsMZ40AGTWjIofWcb63BQCLCHLa/xs5xA0JYk9zhcBRamwZH74MCtmLnCspsHFXgrFgwrgApLpdp527+wyu2iDSttdGjjbaWjBX86kkMjyDHSImMRkQYrJBKyQJ6yAKJYOUgAwlvSf7u4H0eJXKAklJlgG42l9mvfc++Um4mnLKj96LayZGrOLFtSFYldiKg28VKECuACneIRJBuV10QZp4OycWCk9CzLsdfL3C7aSJSPa+lhSzrERWrsc/z7k3PPS6x5LMQYe/bLVuVyGY0vMP7xFO04iMmyQBMsPJF3nHCKhsEKeArqlykkD7+vrw/e//30AwNNPP42LL74478f42Mc+hl/96leRn79r1y4sW7ZM+rtzzjkHa9asyXlNKoALSNp9ltPBLlW6F60MQbY+caADq3IE/Hca7PGKIA08QkDV075ml8QQIKE41juegVsBPBMfeJxn4BlaEDNBPO35TnLBOpauIA1rZhVjTF2QX4IdOo4EVqhdIrupV9g3/tL9LvABPCoMwyjQaiYPKokpIGU6zawA/zarpYogH7g40IH15Gavk7+fk1iTFZckI5a2s37gSU5OCCrk4b3TFG4JJdORatY6aeDuJFfYDFw8PyEwy+Dxgaf7fAczcDPd8qFwI9VY33ZGZiIf57FbgINv5ncx05cD7/mm769vueUWbN++HStXrkQ8HkdNTQ2uvPJKrF+/HieffDJ+/vOfgxCCefPm4dprr8WTTz6JG2+8Ec3Nzfjyl7+M8fFxLFy4EHfddRdqa2txyy234Le//S1isRguuOACfPvb3wYArF69GrfffjsOHjyIW2+9FVdeeaVv+1oeo6OjuOaaa7Bx40Ycd9xxGB0dzcvHogK4AHvSSYg2XKoIYpmOJ9p0hiVwsorf+wHMt+3d1rOeHGFgSTmejfJNqjzHNRwXCtsojmaRxGRrYxp4oQJ4TNPSw6I5Bp4Orrw3XPpanQV6R0LRdRJIHpijhu/ymG94E92lS2a++c1vYv369Vi7di2efvppXHrppdiwYQNmzpyJM844A8899xzOPPNMAEAikcCzzz6Lnp4evP/978cf//hH1NTU4Fvf+hZuv/123HjjjXj44YexefNmEEJcLWYPHDiAZ599Fps3b8Yll1yCK6+80rd9LY8777wT1dXVWLduHdatW4eTTjopL+etArgA9sUoVx940JBhfno6K7kPY+C8M0LWn8PaykfphWKxTCc5GRx83IOArecwCSWTYhybgTMbYYE0cFY2L97wxlKGPcpNC5BQUlxxFZNQwm2EhW26xlorx0KuEQ8CmPJk4dRTT0VHRwcAYOXKldi1a5cdwBk7fvHFF7Fx40acccYZAICJiQmcfvrpqK+vRyKRwPXXX4/3vve9Lj39sssug6ZpWLp0Kbq6ugD4t69dsWKF/brVq1fj05/+NABgxYoVrt/lAhXABVgTZjRXP+dyghHkA+e+iE7TLhbggqUMpzBFHOgQrR+4aLMLK+ThmTqlbg08MwbujIxLGZkNRM4E7PxSgoTCGocFMf+Y4EKJMxdKoI2QFfIULonJ+qfzMlS5oLKy0v5Z13WkUin73zU1NQCstrPnn38+7r//fs/rX375ZTz11FN44IEHcMcdd+BPf/qT531Z29qobXYznTgfBSqJKYC5KsqXgUccqZYelmDb7PwklHSlpTjXEnCqJSP3AzfchTxB8xz56kXHB565jZBv4FXIQh6x3J9vFRCW6I2nG6eluHNmHnI/GNx1WiiSwc/dZMcsVdTV1WFwcDCj15x22ml47rnnsG3bNgDAyMgItmzZgqGhIfT39+Oiiy7Cd7/7Xaxduzbwffza14rPue+++wAA69evx7p16zJaqx8UAxdgmO5CnnIL4MFJTHdfEF4qCirkYc8R+26kOMYYBhaQ+EKeIAbufm/rb+Bo4Bk0s+LOr1DNrKw1ETsRqaXlEpb0DdPedY6pA46DJXjYBeUCfeFcKO5rpHS/Cy0tLTjjjDOwbNkyVFVVob29PfQ1bW1tuPvuu3HVVVdhfHwcAPC1r30NdXV1uPTSSzE2NgZKKf7jP/4j8H382tfu2rXLfs4nPvEJXHPNNVixYgVWrlzpCfDZIqcATghpBPBjAMsAUADXUkpfyMO6igLW7S7GuVDKMYnpx/ZYDDHSgSaWDqLsdTIwFgbAw5jZFzoSA7d7obDXBA8tkCUxbQ08IxeKc9NKGtHK/rMBa9bFf/5OB0YzsH9LLN1R0Hb8pOWXIMmCtTBgiVJKad636OI1UsoMHAB+8YtfSB+/44477J/5oAoA5557Ll555RXPa15++WXPY3fffbfr32Hta+fNm2fP2qyqqsIDDzwQeg6ZIlcG/p8AHqeUXkkIqQBQnYc1FQ18HwvbsVHCrEOGoIEOTAZJ2rqsxrFyv0Ie5/1iwgT0sDapPFhA4rsRxgNkApn1jtkIs5+JWTgbYUzXMDJhpHMo1jHimuUkCdPA2fxVV4I3xB6YNNw9e5IGRUUszwGcOV3KUAM/WpB1ACeE1AM4G8DHAIBSOgFgIj/LKg74tp/lXEofZbtub4/DGDhnExRL33mHShjEfuB2kPK1EXolFFaJmQmL5ht4MUZZCDA5iB+oYVsnjWANnAVr3gfOj5GTwdHAC6dPM695OWjgRyty2U8uANAN4C5CyOuEkB8TQmrEJxFCbiCErCGErOnu7s7hcIUHX+Yd5o8uVYTpvGy7zipONY1AI/5SEe9qERsnMfacWT9wTiYIkFDE3tgAV4mZEQPnJvIUMolpSyhOz3E2Mi5slFs87fl2Ps/gqfSUOgOa7WZdBdgp8tOpgPLbjR4NyCWAxwCcBOBOSumJAIYB3CI+iVL6Q0rpKkrpqra2thwOV3i4Pc/lp4E7X+wgyxobPOBY6oLmLyZN/+dlxsDd/cDDnBbuqsTsbYTO35FN5CmQBs6dn/uGZ6aZbPBNlVJggmuXy4K/zKKWMifnOrVZfplo4EcjcrmaOwF0UkpfSv/7V7ACetnCSZyFVyiWItj3K4gR835sFlSCepJ4nscFXL57Xhj0dPBPcTpvkIPC4Niow8CdTn9RwbtskhE969mAySBJg9pJVrbDiHJTBbgbVEjQdHUt1INzGLkgmb7h6WVIZo4WZB3AKaUHAewlhBybfug8ABvzsqoiQayi4x8rB0RxhbDudby7JKitKx98xG192KQZHiyojaecako9zTxNybH5hLKeEwO3znHCMEEjetazAUtYsuIXgLWJDdfAPe1yucSh7O/i9sgX7jo1TBNxzlJbTt+FowW5ulA+BeC+tANlB4Brcl9S8cAPlY0XkNkUClEkDb7tKXte0FQXXmqJC9WBYZNmXMfV3UGK77GRNE1Uarr7uJylzltKn0EAF45b0JmY6X7gcdeN0ZJQAm+qosbP9aOXBU0xEQwUhh0nDXcvlKAB1ArFQU4BnFK6FsCq/Cyl+HAlMW3mVj6sI6zrHcD1nuZ02aCEmVjIwwf6JKdTh4F9nqO2lzu4xwYLSNa5ZN/Mih0jm2k+mcCea8lp4HHNkVCCbjrsd/z5scAsc+m4rJgF7JppmBRVcT3UqaRQPKhKTA5ir2qgvBh4WN9p9rtUenwXY6NBDZFYxR+AdHMmCQOPEBQ9MogeLBPwlZiEMAaeOYsmxDoOm6dZSBcK08BjnAuFJSeDbhzOZ8O7oPydH+LkHqAwDDxlmNArYxlr4N96+VvYfGRzXteypHkJ/unUf/L9/T//8z+jtbUVN910EwDgi1/8Itrb2+0GUlMVqhcKB9HiBpQX64hi62O9p3mmqAeUbfNDi8XSd37wcBjigpTB2wNlgYGvxBRthJmyaF0jWVkQMz1GKl1NyecMAGvdgb1QNPfEoRifOAzQwAvd6pWf+uO3llLBddddh3vuuQcAYJomHnjgAXzkIx8p8qoKD8XAOaS4xBn7UpWTCyUKI7b7bhiOIyOe9obLwA9s0IVSeud40VwoAB+keAbuPbadxOQ08Gym0lvPJxzzL+REHupqWcvftBqrdd/X2jr9hPPZBMkWTtfCwtoIU1lq4EFMuVCYN28eWlpa8Prrr6OrqwsnnngiWlpaJn0dkw0VwDnwGnJYgUspwq0byxHXnfFdvAsliIFXxnkGzmvg0Rm4xyrHOygkn7EtZ0l6oQT1FZHBYuCOw6MQ4G2SMe6GB1g3LV2rCFwfe55GWCMs/yS64zbibIQFYeBWJaZOymM3ev311+Puu+/GwYMHce211xZ7OZMCJaFwcAbpavb/y4mBhw3PBZxgzfdMYW1PZWBeYPa+/Jc4Ew3cSSY67WntRJ1MJuAqMWNCki9TGSSua9xA5ALZCHWuaIe5dnT3OQetD7CYuqOfR8gPaMS+IRWEgaedSuVCZi6//HI8/vjjeOWVV3DhhRcWezmTAsXAOYjNmeIBga0UwXf684M9f5FzRojMWnxPO4mpCaX0rqk5weCtcjHuuIBPkOJ6oThJzOwCuEsDL1ghD4FJLddSIs7Ozznn4L9J+gY1YXCftf/Njd9pOQy8MBJKlGrdUkFFRQXe+c53orGxEbruL1lNJagAziFpuJOAYuVhqSNKd8CYTjCRYo2d5Mza9Z4uF4q7Q16m3QgBSyaIawLLlMoE1ntrhA+EmXcjBNIMPOkkqAsB3q4Yq467HhtLGoE2QlteShkQE6CyHSDv0Al6Xq7gx/OJnShLEaZp4sUXX8SDDz5Y7KVMGpSEwkGUBIImxpQiMinkSQUwa9d78qX0gg88lYkGzlnl9EgM3EqyMhsga3XN5JdMMBkuFHYzHE+ZLsnJeSy4OhawPps45/gB/Bg479ApnITirQEo3e/Cxo0bsWjRIpx33nlYvHhxsZczaVAMnAPfjxlwdM1yQZTKSNY1z6RwbY8DbYR8ZaGRJQPnkpjO+wXIBKa7/DyuaZgI6HUedmzbB17AgQ6Am23z8kbUXih873VA3s7YPfy4cBKKqwYg4CZfCli6dCl27NhR7GVMOhQD55Ayp0YSM1Bv1TWPHhxLtzOVgQ0OAFjJvaSdbJSp9JweHPMEKZmN0N3Bz7mpZn7JxiaDgfOeb4FFA+HefPZah4H739z44cdBBT+5wlUDoIdr4FGH+yr4I9PPUAVwDnwlJsACW/lclPxQXD/EOEsdz66iltK7NPAItkX+uABzWrglFL8g5QqAgjacCWKaZjPwgs3E5J0kmvv8xJ896+N3J7r75ib1yHNNywrZKZC/eQclugEgkUjg8OHDKojnAEopDh8+jEQiEfk1SkLh4C4eYc2byoeBizcgGayycvdkG7FAh4fV0MhhlDxbzmiosa0RG6itjLnWKa82dM+vjNm7ouwklEK7UPgkrbjDAEK8+VyfmKaaCtf7ST3y3HVayLbHrmrdEA28o6MDnZ2dKPWhLaWORCKBjo6OyM9XAZyDmJSL6aWt+4mIokm7HRmO1h9kI4xzn4fMhRKF1fL9PhqrHYkK8C+lj+eNgRPOQ15YDZxS7vrhNfyAmxz/2TgOnaD8gLeUPt/DFuzhIJwcFPRdiMfjmD9/fl7XoBAOJaFwEAMg6zBXLogyYEHmiRaTkzxYObX1PPeUGKfwKTyoxl1JTHcwlskEKa6vtvV6tzacCdjEG6CADJzX6wUXChBeHSu+jz0qLdBG6CQx832dehxZevCQZYXiQAVwDinBB16hEyRT5SOhRPKBc1thWyoKSFBZBT9ya1sUyYaBl0uiaOApk7p6nojacCZwB9dCSSjOWp2Ratz6AzVw7/P0oM/G5ULx99LnArGqN2j8nULxoAI4B+9F6z/ooBRhmOGJOhkrZJ30ZOCHEYjl3UmuWjIMcVeQcrNpmf6eElwouhDYMgH/moL1A+feNy652QT7wL1ulaCBDnzuwSmGyi87FslAUL8cheJBBXAOIjON6f59sksRyQiFPDEJq2Wd9GSwmLCbMadsBh7d2cE/R0zySXuhGF4fOP//TMDfPAo5kcf52bvWQBuhZIdgT9qRyUtcO1k7j5Dn4Oo4mvgEdvl8F44WqADOgf9iAMGjxkoRUSQNN1PkCnR8rHzWHEmh8MZwM/BoMzElckiABm6Ych94qTJw903Cu9ZMb6pBLhR+p+jYCAssoSgNvCShAjgHfpAuUPrVZyKcbW9AElOyrWeDjr3v53bliEUjzGYWpbRdd7FMt4TiF6SCpIVMwL8mGwYfBbKbRDyihBKXfTZBEgo/eIQlO/POwN3XkjhOT6E0oAI4B6b3soAUj1nl2+UCQwi4Mki12gAGzr9G9G0nzeA2qTxclkDGUAMlFNMlO8SFNWQCvjI1017iUcF/5vzOxnnM/7iy50UZNxdL94qJBeQwsoVTZVueZOZogQrgHFJc4QKQHkpbRhdtNBeKJJno40JJCpWWHg3cCO7xwUMmJ8Rt9ijTed3vzfcuzxQyP3m+IXOh8LJKUHsD9/PYZ8PkqgAbIcfW8y1vOK2CuUS3klBKDjkHcEKITgh5nRDySD4WVEykDCowxalXiSkr7/ZjcHzPDet5bg2cL/QIQ0ziQgm0ygnvLVoZM0FULToXBCVpw44rk4r0CBIKL4HlO8EoFmnFlQZeksgHA78JwKY8vE/RkRK37WXWCyVpZKaB8/3ATQqYwrl6AoXuZsypDCQUd5By+6TlGrgplV2yYeCym0e+EZcEa12yfhnEXZ/1f393ibfgLP/6tNjYTTHw0kROVzMhpAPAewH8OD/LKS74KTWAf3KvVGHb+gKCRbAbRAjgXMUf4GXM1lCIiAFcEuCCGjaJ783b2TKF7OaRb+gSFu36rAPWzXRs633cuxNpctmg9uxM673zPy1HbFRWCJ1dIXfkSke+C+BzYGPDJSCE3EAIWUMIWVPqjW74snGg/HzgUTRw+VZf3nfD+RLLXSNiz+4guK1yQpDyS2Jq7t0Qv+ZMELWgJhe4C5W8EkqQBs4/VzxPPwYu7hQLZSPkb95KQik9ZB3ACSEXAzhEKX016HmU0h9SSldRSle1tbVle7hJQdJ0B42KMusHbkQp5AlgtWIy0dlGi0ElLaEI1ZJBEKfq8P83JDdJvo0te724/qiI2tY1F8gYeNRSev657PyCBgmnDLd0FdTLJluI81WD2i0oFA+5MPAzAFxCCNkF4AEA5xJCfp6XVRUJhiihBFQoliKidAcMKhoRA6n4fmJiLSm4dsIgMlN7nmOUSkyu9W2miElcHvmGNDksKUTyfb2kg6GfNCJ65OMFaLrmtFZWDLyUkXUAp5R+nlLaQSmdB+DDAP5EKf3rvK2sCPBIKOlRY+XSpD5KYY10SALrKS0ycCEpalvbOBthJoUxYgEP035lrWxTXBtbfq3xLAIwz/gznacZ+Ri6/41R/Fn6eonG7zdIOGW6k+1hwxaygWe4iUbKajd6tED5wDkkDROykuhy0cGTZvjMSF1gboATFD0auFDMIU5AT0U4Hg/ZpJqgMn6ZLJHtRB7rPQoTvK1jeIM1G8gMRNDAJT53XZPnYDzTivTC2Qj5z10x8NJDXgI4pfRpSunF+XivYiIl6K7s4i0X5hGlsCaow5+oo4r9MNgNzXahCJJTGMQKTOs9Nal+mzSoT2IwGwnFK0/kG+LOzfk5mnbvx9r9G325n5f3JKYwnk9p4MH43Rv78fy2nkk/rprIw8GqxPTqruVSjSlWksog70Yo9xw7pfSitc2xEWbCiGVd+vyYnZiokwW4qMhlmk9UuM5JkH7GEb7uuO7dJcR8mqnxw4bZa/JuIxRu3oqBB+P2P2xBW20l3r6odVKPqyQUDuIYL1tCKZMmPobpZq0yuItGgj3HSUnFH+B8uUUtNgyyPtl+gwKSfhN5spBBZM2l8g3XWgMsk76v19w3U4Axaz+Lpft4ee8HLuQ/Sk0Df2NvH8657c84PDRe7KUAAAbHkth0YGDS82UqgHPwKx6ZSgzc1SHP02VQzsDFwGvYNsLovVD414vMWjq4V5zII9HPo0LWXCrfkPVtsY7tDcxBrxfrEGSsV7RYFrISs1QZ+B82dmHX4RGs6+wv9lIAAANjKQyOp9DZOzqpx1UBnIPlfPBKDKXEPIJgRChtlzNw+Y1KVo0HOEndTHqhWK93mmfxj4k3DkqpJEh5XxsVshFn+YasyyN/7PAbq8SFosntgaLF0i/ZmQscBu7kHvh5qMXG2r19AIBth4aKuxAA4ykDE+nRixsPDEzqsVUA5+DpRlhmATyKJi31gfuUtPOjuwBJKX0GvVCsY8sZuMjsZKPacmPg0RKJucBvaETUFgBiPoK9j7+NUEgEF8hGGBNuLKVAwk2T4o3OPgDA1kODxV0MgMGxlP3zxv0qgBcNfsUj5ZJ9Fws8ZJD2JLErLIObWYmfR6YSisjk2WPiDTJluo/LrzuXfuCFlFD4fiayXVwY+5dKKFFthD5aeS5IhlhIi4kdPcN20CwFBs4H8E2KgRcPhqd4pMxshJFcKM4XkhW1+Gn9fs2sUrYPPDxpysNm8p4iFB/7okTuyWampa1DF1BCAeRySWQboX1+4Rq4t09M/i1+4jAPu1q3BMgMk09OW9CMrYeGii7rDI4lAQA1FbqSUIoJUYKwk3tlk8Q0I7tQZBY9r4Ti7ofhGeiQYSm9M6xA0Hl9tHcZk81lJmYhGTjAFUaFTOeRIZNSelmfmHz7wMXWxLES2o2u3duLusoYLjx+OgbHUugeLK4ThTHwk+c1o7N3FP2jyUk7tgrgHJKeJCa7aKcOA5cny/wkFG8iix0HYGPPogdF6bBfSbm4KN241pBDJWYhbYQAf6Pw9l4JL6X3sVhKAnNyErpmiq2JS42Br5jdgGPa6wAAW4sso7AA/rb5zQCAzZPIwlUA5yBquowpTqSKf9FGQRQNXDaaLExC8TSfMpyhxtlo4GKiLqwPOb/eUi3k4Y8jm7AT5p6RdTAMajPgavmg5X9aTlK4eQf1J59MjCUNbD4wiJWzG7F4Wi2A/OjglFI8uGYvugbGMn4tk1BOW2AFcFEH39I1iHO/8zRe3HE453WKUAGcg1hKH4+VFwOP4kJxfMleCcXDhIVKTJGFJY3MbISyYhxZEYpYRMIfu1RL6V3H0b0kILwSU/Jan5F+ScP0MPCCuVBCEt2TjfX7+pEyKVbObkJbXSXqErG8OFG2HRrCZ3+1Dv/+aObDxRgDX9hWi5aaCo8Ovr9vFDu6hwuyA1QBnENKSA7ZF20ZaeBhQUrmyBD93QyG0A/D0cqZBp6ZjdDPheKxEQpFJEBuMzEno5mV6zgSCSW8EtN7fn7FM562xwXxgYvj9NzyWaY4NDCGLV25B1qWwFw5uxGEECyeVitl4J9/6E3c/Mu1kd/3qc2HAAC/X3cAB/u9LHw8ZfiunwXw2soYjptRj00H3M87lNbop9UlIq8nKlQA5yBWYrKgMRVdKLIJOX5uEF1gsHwvlExYrSzAydijlIHruTPwgksokkpTp+Nj8Lrjkh1G3KfHidizJ1aAgcNMjnOcSsEMfGAsiW8+thlD4ynp77/95Fv4+L1rcl7X63v7MKuxCm11lQCART4BfF1nH97c1x/5ff+06RBmNVbBpBT3vLDL8/sH13Tivf/1DAbGvAnKwbEkqit0xHQNS2fW462uQVfMYElWtuZ8QgVwDqItTmScpQ5RApJBts13GLhPIU86WLC4ZJfSRzie69iMjXp6ofi0sQ3wrGcCmT+7EJBNDRKHYfjBLuQRGHjUgc/5Jhni3zZMA//x6h34wV+244Xtcp330OA4DvSPZW35G50w8N0/bsEfNnZh1bwm+/HF0+rQMzSB3uEJ1/OPDE+gb2RCfBspeocnsGb3EVxx0ixcsHQ6fvHSHoxMuG9E+/pGkTQoeiSOl8GxFOoSVl/A42bUYSJlYmfPsP37QwNjqE/EkIjrkc83Ko6qAP69P2/DQ691+v5e7G/NvvCsTLbUESWpGGQj9M7EdG+jCSGI68Qu7864ElPix9YDNHCZJz+rkWqT0AsFkO9uos7y1AW7HvvZb+Cz2BKhEEONZechO87AWBJ3Pb8LANDj01yqdySJiZTpy9D9sL17CN/78za86/a/4Lt/3Irzl7bjixcdZ/9+EUtkdjssnFKKw8MT6BtJRrph/GVLN0wKnHtcO647az76R5P49Wv73OtP3yB6RyQMfDyJ2korgM9rqQEA7D0yYv++a2Ac0+rzL58AR1kAv+u5nfjVq/IATim1elBzX4zWGmvLc2gw88x0MWB9scMq/rwFMX6jzWQFNbwuG8X14j62xEYoLeRh2rusojF7Bp7NazM7jvcmE7WUXmqx9HGhiD1o2N8knwUtYrk+u65kUs3PXtht68AyhgoA/Wk2fGTYy4pNk+L6e17BM1vdQ8+v/unLOO87f8FtT7yFGQ0JPHDDafjeX53kCoaLJE6UofEUJlImUibFYIQbxh83daG1thIrZjVg1dwmrOhowN3P7XQ953B63TJWbzHwOABgZmMVACtxyXBocAzTCiCfAGUcwPceGcHj6w9Efv7QeAo9QxM4IElQAE6PB/6L0VAdR30ihj3c3ZTHyETKdacV8dc/fgk/f3G3Z92HsrAqRUEURizTg+1RaaKEYlBoxBqwa78+zZitgJGZJs2eKzpg/Ee5eXdD2c3EzF4/z+o4/A0vov4u86r7SiiGKRQL5V/qE2/OogY+ljQAWN+BHz+zA+88tg31iVggAwecQMjjyMgE/rjpEJ7d6gxESBkmVm/txnuXz8ALnz8Xv/rE23HaghbPa2c1VqEqrmNrlxPA+ZtE33BwUU3SMPGXLd04d0kbtLTmf+Hx07G9e9glozAGLrsB8RJKW20l4jrBfi7OHBocVwGcR/fgOK760Yv45H2vYTxlRHrNnsNWoN3fNyplKmLva4a5LTXYc0TeIvJ//rIDF//3szB9yp2f296DV3f3uh6/7p5X8M5vP43/fWVv3kuAU0IPbRmko7t8AoCs2yBLOsr6lYQhOgN3N1LiX5PLVPrCM3DimUkalzwmfa3E6ujn7/bMbi3A4BHRkcVr4K/u7sWSf34c537naXz83jXoHUnixnMXo7WuEt2SAG6Y1E7+HR7yBkCW5Ovhftc7kgSllrd6RkOV7zo1jWB+aw129DgB3P0+wTr4ml29GBxL4dwl7fZjLNnYM+i89ojNwOVJzPo0A9c0gvb6BA6kGTil1ArgSkKxMDph4Pp716CzdxQmBbr6o5XR7jliJRXGU6ZUx5IVjwDAnOZqX5a9s2cY/aNJ9Ax719AzNA5KvZpgZ+8oDErxuV+vw9/+7NW8ugcy0cDjEgYua2Ylvh/b1st06jDIOu7FJCPV2M2UD9bzWmowq7EKC1prIh/PWXP2FsRMjyPe0GK69zH5ayUJUB9/tyfZbktg+cvViJ05+YT+5oOWz7mtthJrdvXinGPbcPLcJrTWVrqCHsPgmBWMAeCI5LvCbHaHud+xgNlSG85cZzYm0DXgfS0QHMB39gzj1ic2o0LXcNZiZ5JOW/qY/M3oCJOAfCUUZ7jZzMYq7O+zGPjAqCXnFIqBl91ItX94cC3WdfbhY2+fh7uf34X9/aOY01Id+rrdh50gvL9vFM01Fa7fGxLrGgDMbq7GkxsPSi16rGqrs3fU4/FkFxTPOEYmUhiZMPDZC4/F0HgKdz69Hev39eOE2Y2h64+CKLY+cQoP/7PoMJBp3DFNg2FQx2KYjY0wZCIP255XcEFqekMCz91ybuRj8ciFvWd6nLjkhic+JoOs0jTuq4Gbnl0MkG8GLnrNHQ380MA4CAHuu/5tMCiFnt5dtNVVYpOknSpPmHoCGDgfeNmkHfF7KsO0+gRe29PneS0gZ8wA8KPVO3Dbk2+hMqbh1itXoKbSCYWt6QDOyFfKMO338dfAuQDekMCa9M6b5c8KYSEEyoyBdw+O49E3D+KT5yzER0+fCwA40B9tAgavY8t0cFnxCGAx8KRBpcdhzEE2hYMFd55VsGDeVleJ85daWzaZppYtojBwS+eT99qW2flECUXXLM2aBfuMGHjEQp43OvsR1wkWttVGfu8gTMZEHsC6drwM3PuY9LUSq6PMoQOkGbhMQskjAxcJi85dI4cGx9BSU4mYrqEyptvHb6uVSyh80JNd7yyA82SnJ/281trwAD69PoEjwxO2nMrr7LLjdfaO4OuPbsLbF7bgqc+8A5edOMv1+9Y665gsgPdxzal6BU09aZgYTRp2EhMAZjRWoWtgzLrZFbCIByizAM6C8Kq5zZjZwLK90RKCe46MYEaD9SHKgrGseAQA5qbZvZjIpJTaQXqfJICzROXhoQlb62YXd1ttJZqrrYsknwE8igYOWMGCv1FpGoFGvAxOdkNgLU4Nm4FHD4piQRAgb4X68s4jWDarAVUV+fHNTsZEHuv9icdr/p5lM3DtmfPDX+vT6EsMyix5rMssfnlk4EnDlN/kDROHBuRJudbaCgyOpewdFAPPgmXXO2OpluxonQNj0S014cy1vT7tFuN2vYm4BkLkjJlVWl5zxnypNs2OyeQgfs2ihDLEVWEyzGyssjzjQ+P2uU2rLzEGTgiZTQj5MyFkEyFkAyHkpnwuTIbOXiuIzm6uQlWFjsbquMuuE4Tdh0dw0pwmK0MsCfoy3RWwGDgAjw4+NG7JIfy6eDAJJWVSDIxaf2TGMFpqK9CU3hqGJVkyQVRfdkzTpNKIrBuhnwbOLIeZ9Nj2k2/4wDOWNLCusw+nzmuO/L6hx83BgpgJZHr3GYta8ffvOib0tXEfecnXYimRN/IZwEUGzmvgVlJOFsDd0gND36h1jTdUxaUuFcbAx1MmhtPfqcNDE9A1goaquOf5ItrTQZgFyyPD42itrURDVVya7+rhiJQMFTHNtVYWwOsSMc8NgdknRQkFsKRadlMpRRdKCsA/UEqPA3AagL8jhCzNz7LkYE6SjiYrqM5oqPK1BfJIGib29Y1iXms12usTOChh4OIAX4YZDQnoGvEwcD5pEiShALCTnOyCaK2tRH0ihphGpLaqbGFEaGYFpBm4EHit4cISDVx0oWhWgyXGYuqroqdRfCevcyxz7d4+JA2KU/IYwGUdGAsBGQOPirpEDLpGUMVV67EeJ7xbSWaxtGe35lFCSQp/e14D7xoYQ7tEErDdG4LOzRj4graaQAkFcJj34eEJNNdUuCysfmAB/GC/89qW2ko0VVdICVKU0vbW2gpPAF/YVosjgoTC3DW8hMK84Af6x3BocBxVcd3F0POJrK9oSukBSulr6Z8HAWwCMCv4Vblhb+8I2uoq7ZLUWY2JSAx8f98oDJNibnMNZjZUuTyaDOLwAoaYrmFWY5XHSsgkkqbqOPZJ1tA16NW++cQMIQRNNRWeEuBcELWwRteJZ6fh19ZVfD+mWT+/zfLsnjrf6831g18/cJ45vrLzCAC4yqVzRS6NsDJBVMeJDJedOAu//sTbXYGABVD+zyKzWPLBNV8QG5Wx8xpPGegZCmHgQjFP70gShADzW2t8bYQswDFCc3hoHC0REpiApYEDXN5paAItNRVorI5Lk5jdg+PQSHCCtLW2UhrA+0YmXDdUxsDrXQzcKeZhu5UwG2m2yAslIYTMA3AigJckv7uBELKGELKmu7vb89pMsPfIKGY3OZ7QqAycOVDmtFRjRmNCqoHbg3QlX8C5LdUeBn4wfbGcPLcJnb0jHk/3oYExzG621soCd8/QBOoqnZ4IzdUVGWngA2NJaac0BkPCmGWQSSiywbgpoW2p9Twr0D+7rQfHz6yP5BJg0KXFKu5p5y/vOoJj2+vQWB39fcOPG62YJld85G1zcNN5i7N6bSKuY6XgRpLNoUwJHSIBjoHnsR9K0pAX8nQNjMOkkGrHrXVyCaV/ZAL1iTja6ipxZHjC813pHhzHsdOt4Qw22RmeQEuEBCYANFbHURHT7AB+ZNgK4M1+DHxoAs01lYHXQ2tdpb2TYCRrQVsNUiZ1tQNgP/M33vqqGKordOzvG8OhgcJVYQJ5COCEkFoAvwbw95RSj4eIUvpDSukqSumqtra2nI61t3cEs5sdy+CMxgT6R5OexjMiWPCd21KNGQ1VONg/5im+Eadw85jdXI09h4ddjzEJ5cQ5TRhLmh4ppGtgDEtn1ANwMuo9Q+P2RQ5YDCATDfwrv92Iv/6J5x5pI6oGXhmz3AM8ZG4QmYSiawQDo0m8tqcXZ3Le2SiY0ZBAfSLmSk7y085ThonXdvfilPn5Y9+AdXP67IXH4qLlM/L6viJOW9CCK07uyNv7yabgONep3Eb46u4j+NPmrpyPLY5tYz+zAhVZUGKMWRxx1jeaRGN1HC01FZgw3P1QRicMDI6nsMQO4OP2/6MkMAGrR097fSW6BsbSfVDG0VxbgcbqCl8GHuZuaauttHcSh4cnUJeI2efMO1EGbQnFYeCEkLQXfBTdg+MFc6AAOQZwQkgcVvC+j1L6UH6WJEfKMHGgfwyzm5wAHtWJsufICCpiGtrrEpjZmLAyxJy9j1JqT8uQeYXnNFejdyTpaiXZNTCGusqYPdaJd6KMpwz0jiSxZLoVwB0G7t4WNtdUZKSBr93bi509w9KOcKZJYdJoLPNbV6zADe9Y4Hosrnl7SssLeTSs6+xH0qA4c1FmAfySE2biuVvOdd08ePa46cAghieMvOrfDH/3zkU2yysXyCosWfJYNvXnQP8YrrtnDf7+gbU5s3Fvb3zr533p71q7hIEn4rq0nL53JInG6go7IPMyCgv2S9Jkx5FQojNwAGivs4p5BsdTSBoUrTWVaKqOSwlSz9B4qC+7ra4Sg+OWo6Z3xNLjm6q9xgNZEhOwyMqBfktCKZQHHMjNhUIA/ATAJkrp7flbkhwH+i1fJZMlAATaAnnsPjyM2U1V0DRil+UeSF+Ie4+M4KM/eRlf+/0mnDy3yZ5rx2OuxIlyaHAM0+or0ZGWdPhEJss8z2qqQmN1nNPAJ2ydEACaauKRNfCxpIGdPcMwTCqVjWSNp/xw5uJWj8dal0xAl7WLjaUllIqYlnGg1TTi2mry6zVMipd3Wfr3qZK/wdEIx5/vBGNZwRmT/b72+43oG0liYCyFV9KfJWBdO5nmWsT8ByMG+wMYOOCWHhj6RybQWBVHczog86Sle8i6lmc3VaG20gr+4ymLlbdGqMJkaG9IoGtgzP6uNddYTq+RCcNja+weHPd1oNjnUet4wY+kE6pNNda16w7gFqmrFQL4rMYq7OgextB4qmAWQiA3Bn4GgI8COJcQsjb930V5WpcHLHi6GLik85cMuw+PYG66zaMY9G+8/3Ws3duHr156PB7829NRXeHNFs+WBPCugXG01ycwyw7g7uAOWCylpabCLuYRdb3m6gr0jSYjJZ+2dA3ayay9EtsiO58opccyxDXNNevyJ8/uxCu7jni2seyLvGpuU176GzNdPGVSvLzzMDqaqgJ7XxxNkLX5lRWcMQbe2TuKT5+3GBUxDX/Y6MgoX3j4TZx3+19cFYph8EgoLICnrzM/VtkqKeaxGHjc7u7J5314R0hzjZUTYr/PJL9iMfAxu1S/pdZKYgJuHzqlNBIDZzeP7sF0AK/2Z+AVEklyRkOV3QmxJCUUSumzlFJCKV1BKV2Z/u/RfC6OB9OxeQ28vT4BQoIlFEop9hwZsf3cTtAfQ2fvCN7Y24dPnbsIHz19nq9laY6kmKdrYAzT6xOoT1gdC3knCtPH2+sr0VJbicNDE0gZJnpH3Ay8uaYClAL9o8Ed0wD3oFSZbXFdZz8AYPmshtD3kkHXCPYeGcH3n96GK+58Hl99ZCNOX9CCr122zPU89kU+I0P5xA8sofn4+gP4w8YuvPPYaXl536kAp8eJzEbobTR1QkcDPn3uIpy5qBV/3NRlNVIaGMNv1+7HkeEJfPWRjZGPnRQllPTfqW8kiZaaCl+7JK8dM/SNTKCpusJh4FyAP8QF8JbaChwemnDqJTII4NMbKjE8YdiGhZaaSnnAHU9hPGVGDuA9QxMOA2fvx2ngA2MplwOFYUajE7RLOok5WdjbOwJdIzaDBizDfWttpY+rxIRhUvQMTWBkwrArKpuq46iMaTjQP4onNlgs5cLjpwceuz4RR2N13L44rC+G02Gso6naFVRZNry9LoHWWkvnPjIyAUrdpcGsmEfW4EfEpgODqIrr0AjQKWmu9ea+flTENFuTzxT1VXG80dmPWx9/C4eHx/HdD63ETz92CqY3uNkDY3tnZZjA9AMLPl94eD2Wz2rA5y9akpf3nQpgu5MN+/pt54askGdxey3OX9qO73xwJWK6hncd1469R0axpWsI9720BymT4v0nzsL/rd2Pv2yJ5gTzMnAnVAR11msTOhJanQhTaKiK2wH5sMDANWIF3JYay7rHNPRMdpNMk9+Y7sXCM3A+gDPGHybP8I4aFsDrq+LQiFdCEWVBwJJQGAopoZRNM6u9R0YxszHhcUXMbEhINeGL/+tZHBmZsHsIswBOiHUT2N8/hjf29mPJ9DrMi9Dhbm6zYyXsHUliwjDtEt5ZTVXYzblUugbGUaFr6cx7JQ4PHbbLckUGDsBTHAAAT791CMPjBt67wnJObDowgCUz6tDVPyZl4G/s7cNxM+pREcvunvzfV52IgwNjWDytVnpBMrAqteNnZsf0RTCW2dFUhZ987BSphHW0Ysn0OtRWxnDDz17F4mm1+MJFx9nXDK9P1yfi+NHfrLL/fd5x04CHgd+/eQD3v7wH7zy2Dd+4YjnWdvbhiw+/iSdvPjv0c04a8l4oQDCj5MvpE3Hd3l02VseRiOuoqdA9ScyWWsvS11pbgTc6++zfR+mDwmAH8PROtbmmwj42L6H0RCjiARz2v+fICMZTJpprKuzKUFFCEROYAFxEsyQllMnG3t4Rl/7NMKOhyqOBHx4ax1tdg6iu0PHom9bQhwWtta7XbNjXj1d2Hwll3wyzm6uxKx2kbYZtM/Aq7Ot1+owfGhizzfsttRXoHUmiK62Lt0gDuDfB9M3HNuOWh9ZhPGWAUorNBwexZHo9OpqrPRq4aVKrq2FH9kF1ZmMVTprTFBi8AeCT5yzEHX91Yt481UtnNODEOY2455pTM0paHQ1YNqsBL33hPHzriuUwKcWNv3jNvgaD/P7t9Qmc0NGA//nLdnQPjuPqt89DZUzHNy5fjs7eUTy4Rj6VKmWY2H14GGa6103c5UJx/t7tAYyS/Q3FCTZMfmiurXDtOA9xCcWWWksD7+EK3qKCD+C16VoLmYTSPRSNgSfiOuoSMWw5aE2YZ7tlq7rTuSEMjcsDOJNq4zpBU3V4O4BsUTZ0Z++RUZy3xKuPzmhMYPXWblBK7WqnDelt1DcuX45F02qxo2fYxbJnNCbwQto2+O5l0QL4SXOa8Mi6A9hzeIQL4GkG3liF4QkDfSNJNNVUoGtwzL6gWMDe2mVdCDyr8AvgfSMT2Jy+cJ5+qxvLZzWgfzSJpTPqMJ4yPINjd/QMY3jCyFr/zgT5Yt4Myzsa8PAnz8jre04l1FTG8KFT5uCMRa248D9W46uPbAIQ7jY6f2k73ujsx/zWGpy92Kq/eNuCFrTVVdr5EhH//adt+M+ntqIlzV75xmis4ZlJgxkln/yb1Vhld/JrSAex5ppKj4TCJIbmmkoYJsXOnmFUxLSMys/Zd7FvJGnnu2RJzKgMHLD0/C2HrO8hY+Ri9fTgWBJttd6umYm4juaaCiRiWsGqMIEyYeCjE1b5Lm8hZJjVWIWRCcNuGAU4AXzpzHpMq094RjEx//jclmq7gCAM5x1n3Tz+uKnLtgm2cxo44CQXLYeKdYG0pv/wbx20JobwDFzGEABrSghgTYH/7Rv77Qb6S2bUY3ZTNQ4OjLkGLb+5rw8AsKKjMdK5KJQfOpqq8aWLl9rsNCyAX3j8dBACXHOGOzm/dEa9KyHOMDiWxF3P7cTJc5tw9jFtmNGYwAqBEDAdPIiBO9Ns0q1YBQbeWlPhklAODY7ZDJyRm7e6BtGabjcRFdUVMZsJM6dXIq6jukJ3BdzuoXHoGkFjhCZZrbWV2JtuoeEwcHeDrMGxlMdCyDCzMYG2Ak3iYSgLBu50IZRLKIBlb2J3+fX7+9HRVOVbjs0yxO8+fnrki2RuSw0WTavFnzYfsr3i7GJlXvB9fSNY3tGAroExu8iFBey3ugZQoWuujDW7wEQG/vKuI6jQNVy6ciZ+t26/7UM/dnoddvUMg1LLOsl2FW/s7UdVXLcHvCpMTXz4lNl4fP1B/GVLd+hwisXtdXjqM+/AfCG/c9yMejy/fQcmUqYrX3L/y3swMJbCv1y81HfAiK4RwADaghi4UE7P2C8LmM01FTbBMtMmA8bAmWV1a9eQZ91RML0+gcGxIZd7RZQ8WBVmlCZZrC+4tTZHQlm/z7kB+mngAPC5C5dAKyD7BsqEgbPkYYdMA2/0FvNs3D+A42fW+77fkul10Ahw8YqZGa3jvCXT8NLOw9jWPYTmmgrb+8kCOBuEOjiW4iQU6w+/tWsILbVeVtEsaWj18s4jOGF2Az6wajbGkiZ+9sJudDRVoT4Rt29ifCLzzX39WDarvuC9PhSKC0IIbrtyBa49Yz6WzgiXsha01Xqut6Uz65E0KLZ3OzMkx1MGfvzMTpyxqCVwOhRj/UEMnAU6FsBZ8GRyRkttJQ4PW32/e0cmYJjUpYEDlq6cSRUmA/vO8dp5o1CN2SMU0wXBXXTHSSjp9zPSfVH88kZnH9OWcbuJTFEWAZwV0MyRMHAmh7AS36HxFHb2DAdqtSfPbcarXzofyzNM+p27ZBqSBsUTGw66MvENVXGsnN2IHz+zw7YxORKK9f/xlCm9cMRy+pGJFNbv68ep85uxam4TZjQkMDiewnHpUmN2s2CJzJRhYsP+fiyf1ZjRuSiUJ6bVJ/Av71ua9bCLpTMsyXAjN/rs4df24dDgOD7xjkWBr2WsP8hGyJJ/zK7XPzIBQmAP/W2pqUDSoBgcT3EecDfZATJLYDKIeScAnpay3RmUtrMbS1wnqEvr8Y3VcYynTIxOGHZPF5kPfLJQFhLK3t5RVMV1qa2ora4SFbqGbekkIdP3ghg44NxRM8HJc5vQUBVH/2jS1QuCEIJvf2AFLvqvZ3Hz/64F4FxM9VUxewiCjFWIF9jre/qQMq1+2JpG8L4TZuKHq3fguLRWP73e6k/OZKWth4YwljSxIgcHytGK9T3rsWtgV9avX9ayDPMa5uVtPZOBeS01qIxp9vfENCn+Z/UOLJ/VgDMWBbcGZr70sDL0Nq6cvm80iYaquC1Z2In7oQk7yDMJpam6AtDGEKt9C0OxfXhkh9wt44fh+AHE6rvQbXbhkR3bsKxlGRqFds89Q+OR815MDmqqdnbO9iQtrq2sn4QyGSiLAH7q/GbUJWJSvVrXCN61dBp+88Z+fP6i47Bhn5Vhz7dbArCsW+cc24bfrN1v9yBmWDStDp+94Fh8/VHLJcAYOLMSdg3Iu6u11FS4trMv7TwCjVg3CwC4/MRZ+MmzO3Fyuu9ITNcwszFhJ1feTDsKVADPDCY1cf2T12M4ORz+ZB8sbVmKX178yzyuqvCI6RqWTK+z/dKv7bEapN3+wRNC80ExjaC5piK01qC1ttIOzr0jSTuBCTgs+/DwuFNGbzNdDfXtz4M2PokXBoEXnsn8/KpmAU90W/8tb12OhdVfsAkSK6NvjcjA2Y7ZLcmwaswJW7IMs94WEmURwC88fnqgX/sjb5uLR988iMfXH8SG/QNoqakI1OlywXnHteM3a/dL3//aM+fjiQ0HsWZ3r2ub2VJTia6BcVdShEG0Jb2y8wiWzqy3L4rjZtRjzRff5doxdDRW2wz8qc1daKiKY15L5kmfoxn7h/ZjODmMm066CefPPT/j19+1/i48suMRGKYBXcvP7M7JwnEz6vHEhoOglOI3a/cjEddwQYR6CF0jmFYd/r06bnod7n9lL/pHkugbmXCNRWMkZtuhIezotm6evKRRUdWF0YkWfPr4b0W2+DKs3tKNf/nNBtx65QqsH/4/PLbzMZzcahUTGSbF4FgSSYOG7iAY2I6fD+Ds596RCTsHphh4jjh9QQvmtVTjvpd2Y3jcwNKZ9QXzXr7jmDbMbEhILXu6RnDnX5+MV3f32pof4LCOVgkDb66pwHC6Y5pGCF7b04uPvG2u6zmi3DO7uQpPv9WNjfsH8MSGLnz6vMWRsuoKDrb3bQcArGpfhbn1c0Oe7cUJbSfg11t/jX1D+zCnfk6+l1dQLJ1Zjwde2YvO3lE8+uYBnHdceyTPdVwngfo3w5Unz8Y9L+zGb9ftt3qncNIhC9b/9Os3AVi2vBru2DR+EMbQdBzTMh9z6zPri3PpspnYti+Bi5YsgbZ7Bx7a+hC0ij5QCgyMJp2Rhjkw8Ca7PD+J2krLyluocWlRMCUCuKYR/NXb5uDfH90MjQA3nL2wYMdqqIrj+c+f5/v7trpKD3NgF4KUgXNe8L1HRjGeMnFqyECDjqZqHBocx7ce34y6RAzXRZh6ruDG9n4rgC9oXBDyTDnY67b3bS+7AM4S4j95dicOD0/gkhOiubEuWj4j0k5v2ax6LJlehwfX7EXf6ITL3jq9IYH//PBKDI1b/VF4PXrCmMCE1g1zYmlGjawY6hNxfOVSq/nawkYrBoxrBwBUondkwiPZhIHdbPwklDc7rV1wlMHLhcKUCOCAddf/9hNbMGGYoQnMyQa7GGUaOF+N+eibB1AR00I7/bGCpr9s6cbN7zqmqBdQuWJ733ZMq5qG+orsrpUFDekA3r8d78Q787m0goMFzfte2o26RAznHBttUtbn3h2t0RghBB9cNRv/9shGaAQ4b4n7+rx0pXx07u6B3QBMmOPtWbdFZpjfYJGaIXMfgAXoHUnaZfRtEiIlQyKu430nzLQrWQHHDnnHn7ehe3Acl5wwMyvPer5QFjbCKGiuqcB7llvMt+QCuF1p5h/Ae4asAP7OY9tCkyLMD1+fiOGaM+fld7FHCXb07ciafQNAXUUdplVPw46+HXlc1eSgLhHHnOZqJA2K9yyb7ullnQ9cduIsxHUCkzpBLwxsV6Sl2rNi4DwaKhvQVtWGvqTlZPn9ugMcA49eHfnfV52Idy1tt/8d1zXbJnnTeYvxnx9eWdBS+TBMGQYOAP94wbFYMr2+qHdEGY6bUYe6ypg9/IFHc3rKx5MbDuLQ4DjeG6G4aH5rDWIawd++Y6FLa1eIBpOa2N6/HVcsviKn91nUuMgOOuWGpTPqsefICC45Qc6Gc0VzTQXOX9qOR9886HKhBGF733ZoRMO9H704L8NCFjYuxOGJPfjgqg789Lmd1vBjXUN9VW5h7+Z3HYOZjQm8e1lhZ6xGwZRh4IBVav+JcxYW9Y4owznHTsO6f71AKnU0p2WVh1/fh0RckzbsEtFaW4k//+M5+OQ5hdP6pzIODh/EaGo0JwYOWDLKzv6dMGn+psFPFs5f2o6T5zbh9IXB3u9c8MFVswGEd/5j2N63HR21HThtfmbuEz8sbFyI7f3b8c33L8c/XnCMnVDNNT5ce+b8kgjewBRj4KUMv4umoSoOQoCRCQMXLZ/uysgHQdYXRiEamANlYUNuN8CFjQsxmhrF/qH96KjL3zT6ycAVJ3fgipMLu+Z3HNOGn35sVeTpTbnKWiIWNCzAaGoUXSNduPHcxVjcXoehsVT4C8sIU4qBlyP4zmiZ9mZRyA47+i3dmjkVsgV7PXs/BTcIITh3SXskjT1pJrF7YHfON1Ue7O/DZK4Lj59e8JvWZEMF8BJAU00Fqit0NQ9ykrCtbxtaq1rRUJlb9artROkrTx28lLBnYA9SNJXzTZUHuxlM5b+PklBKAOctmYa4rmXdoEghM+zo25EXpsecDlM5QEwWbFkrjwG8MdGIlkTLlP77qABeAvjie5cWewlHDSil2N6/HZcsvCQv77egcYGSUPKA7f3bQUBs/3a+wBKZUxU5SSiEkHcTQt4ihGwjhNySr0UpKBQKXSNdGE4O501rXdiwENv7ttud6RSyw46+HZhZOxNVMa/VNhcsaFiAHX07puzfJ+sATgjRAXwPwHsALAVwFSFEUUmFkgbbTufL7bCwcSFGUiM4OHwwL+93tGJb37a8yicMCxsXYig5hK6Rrry/dykgFwnlVADbKKU7AIAQ8gCASwFszMfCePz7S/+OB996MN9vq3AUwoTl2c5XsGDv8+6H3g1NeQKyRoqmcFbHWXl/X/vv8+t3g6C49SF3nHcHzpiV3wHeuQTwWQD2cv/uBPA28UmEkBsA3AAAc+Zk1/Tn9Bmnozau5j0q5AcddR1oTjTn5b1Wtq3EZ07+DAYnBvPyfkcrNKLlXBkrw4nTTsTNJ9+MoYmh8CcXGLNq81/1SrLVhgghHwBwIaX0+vS/PwrgVErpp/xes2rVKrpmzZqsjqegoKBwtIIQ8iqldJX4eC57vk4As7l/dwDYn8P7KSgoKChkgFwC+CsAFhNC5hNCKgB8GMBv87MsBQUFBYUwZK2BU0pThJAbATwBQAfwU0rphrytTEFBQUEhEDkV8lBKHwXwaJ7WoqCgoKCQAZTvSUFBQaFMoQK4goKCQplCBXAFBQWFMoUK4AoKCgpliqwLebI6GCHdAHZn+fJWAD15XE4xoc6l9DBVzgNQ51KqyOVc5lJK28QHJzWA5wJCyBpZJVI5Qp1L6WGqnAegzqVUUYhzURKKgoKCQplCBXAFBQWFMkU5BfAfFnsBeYQ6l9LDVDkPQJ1LqSLv51I2GriCgoKCghvlxMAVFBQUFDioAF4kEEKKOx4kj1DnUnqYKuehEIyyCOCEkFj6/1PioiSEVAOo4/5dtuelzqX0MFXOQyEcJR3ACSHvIIT8GsDXCSHzaBkL9sRCDSHkOwDWAriNEHIFAJTbeU2lcwGA9LncjjI+l6n2N5GBEPIeQsjs8GeWNvJ5HiUbwAkh7QD+GcAjAAwA/0YIOa+4q8oOhJBY+ku0BMBSAGcA+DWAfyCEnJV+TlmwJEJIY/pclqL8z2UuIaQKwCKU8bkQQuq46+t4lOl5+IEQchwh5GUAXwLwU0LIdcVeUzYoxHmUbAAHcCKAOKX0LgBfBfAcgMsJIfmZRjsJIIScRQj5BYB/JoTMBXAygBcopd2U0icBPAjg6+zpxVpnVBBCfgTgN4SQWgDLALxSjudCCFlBCPkNgN8DOBbAWSjDvwsh5GxCyP8CuI8Q8i4A5wJYXW7nEQFnA3iVUnoGgG8CuJQQkt/x7pODvJ9HKQfwdQDG09LJKICXAaQAXFjcZYWDEBIjhNwG4LsAHgfQAuAfAdQCOJ89j1L6HwCOIYScQCk1S5UlEULYdVIBa/rSuwFsB/e3KJdzSeNaAGsopcsopWsBdME6JwDlcS6EkNNhEZuHATwG4AMA4gAuYs8ph/MQQQiZRQi5jRByIyFkZvrhYQApQohOKX0KwHoA56R36SUJQkgHIeRfCSEfJ4R0pB8eQp7Po5QD+AiAN+B8sbYB2ApgLiFEL9qqIoBSmoI1au5SSum9sLZMKwD8L4CZhJBzuKf/FMA16deVpFaZ/vLPBVADizl8FMAzAFoIIe/knlry50IIWQhgJqX0q+l/n0ApfRDALIENlfq5LAYwSim9H8DdAJoopf8O6zzO4p5X6udhgxAyB9auSAcwC8B/E0KmAxgFMABgbvqpj8A6/+nFWKcf0nmISkLI5wE8Bev7shLAtwghTbCk4H7k8TxKOYD3A3gdwGmEkDZK6SCAelgXqlEGbOJZSmknISRBKe0DkIR1Yd4BK6Az7ADQWQbncxhWJ7VdAPpgSVzfAfAZ7jnlcC47AJxACLkqLaN8nxDyjwD2AviU8LxSPpenABxLCPkPWGudTwi5FsALAP6ee15Jnwch5Fjun/MAPE0p/Qyl9PMAOgH8K4A/wAroSwGAUvo8gHYAc9LvUfQ4RghZQi2MwyKf76KUfhbArbBuPksAvAorWOftPIp+4n5Is4UnYH0Y/5p+uArWdqrk2QSldIz9P836mgB0pbe1KULIvxBCLgZwHYCdpX4+AN4FYB+ldD2ADQB+AeB0ADXldC7ptd0HS9L6NoB3ADgE4AiAOYSQWwgh70OJnwuldB8suaQdwFUALgDQAes82gghXyjl8yCEHE8I+QuA9cwtA4uZzuWe9mUAHwJAAbwG4O2EkBPTv9uEdOCjlJqTs2ovuPN4kzuP+2HdfACL9CwBcIhSuhXATgBvy9d5lGwABwBK6WFYwbuJEPI6gHNgBY5yw3sB/JEFdQCfA7APwGcB3EUp/XXRVhYd3QCWE0IeAvBJWEmxX8LS9HthBcRyOZc/A1gO60uVAvAiLKZ6K6wA+A8oj3PZAyuAr6GU9gJ4BcA4rB1eqZ/HOCwJ4W8B/L/0Yw8AOCNNeJDeuf4OwM0AvgdrV/5fhJA7Yd28Hp/kNcvAn8cNAEApPcTdMFvTzzHS/74H1g42L+dRFr1QCCFxAG2U0v3FXksmSCcrDELINwG8CWv3cB2Ar1NKXyzu6jJDmjF8HdbW/UcA3g/gEkrp+4u6sCxBCPkBgBSl9EZCyE0ATqOUXlXsdWUCQkgLrETmFkrpd9Pa62JK6bVFXlokEELqYRGBJwB8mVL6RPrvUkcp/Uj6OZcAuATADelczPsBLABwL6X0ULHWzkM4jy9QSv9ECKmklI4TQv4GwPsppZcJr7kCwHzkeB5lEcDLGek/biesbfprAH5IKf1jcVeVOQghhN+GE0IWAKiilG4Qf1cOIIQ0wto1nAVgEMC/UUpfLqdzSevaFwD4IqzKywMAvkgpfb3MzuMzAM6hlF5CCJkGYDWAmymljxFC/hXWTun7RV1kBKTP4yxK6eWEEC19w7kVFkMfAPB3AO6jlD6dt2OWyd+4bEEIqYMlmfw6bVkraxCrKClV7HXkC4SQaaXC5LIFIWQ+AI1Sur3Ya8kGhJAZAB4FcD2l9FVCyEcAnAJLMjUA3EQpfbaIS4wE7jw+Ril9I00SXoRlvz0A4EeU0rvzekwVwBUUFIoFtlMghNwMy3L3C1hJy9UATqGUPlPM9UWFcB4rYCXKZwO4GsA3KKVPFOK4sUK8qYKCgkIUpIMegeWg+SiAhQD+XzrhXxbBG/Ccx9WwnCcfpVYlecGgAriCgkKx8T4AMwGcTCl9vdiLyQGTfh5KQlFQUCgqyinhGoRinIcK4AoKCgplipIu5FFQUFBQ8IcK4AoKCgplChXAFRQUFMoUKoArKCgolClUAFdQUFAoU6gArqCgoFCm+P9HrAQF5oXcowAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "THRESHOLD = 12.2 # Anomaly score threshold for an instance to be considered as anomaly \n",
    "\n",
    "#TIME_STEPS = dataset.window_length\n",
    "test_score_df = pd.DataFrame(index=range(test_dataset.data_len))\n",
    "test_score_df['loss'] = [loss.item()/test_dataset.window_length for loss in loss_list]\n",
    "test_score_df['y'] = test_dataset.y\n",
    "test_score_df['threshold'] = THRESHOLD\n",
    "test_score_df['anomaly'] = test_score_df.loss > test_score_df.threshold\n",
    "test_score_df['t'] = [x[59].item() for x in test_dataset.x]\n",
    "\n",
    "plt.plot(test_score_df.index, test_score_df.loss, label='loss')\n",
    "plt.plot(test_score_df.index, test_score_df.threshold, label='threshold')\n",
    "plt.plot(test_score_df.index, test_score_df.y, label='y')\n",
    "plt.xticks(rotation=25)\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\User\\anaconda3\\lib\\site-packages\\seaborn\\_decorators.py:36: FutureWarning: Pass the following variables as keyword args: x, y. From version 0.12, the only valid positional argument will be `data`, and passing other arguments without an explicit keyword will result in an error or misinterpretation.\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD/CAYAAAD4xAEfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABcf0lEQVR4nO29eZhkVX3w/zm1d1Xvy0zP2jPDrAzMMMwIgyyCoBB/YNBEghoV0ZeoMSpJ3hBfXuMWE5cYE+KbGIloFggiYtQhoICyGRiYGWBg9n3vfauufTm/P+49t25t3dXdVV29nM/zzDPdVbfuPaf63vM9311IKdFoNBrN3MNR7QFoNBqNpjpoAaDRaDRzFC0ANBqNZo6iBYBGo9HMUbQA0Gg0mjmKq9oDGA+tra1y2bJl1R6GRqPRzCh27tzZK6Vsy319RgmAZcuWsWPHjmoPQ6PRaGYUQogThV7XJiCNRqOZo2gBoNFoNHMULQA0Go1mjjKjfAAajWZukEgkOH36NNFotNpDmVH4fD4WL16M2+0u6XgtADSaMiKTSdKRCI6aGoRLP14T5fTp09TV1bFs2TKEENUezoxASklfXx+nT59m+fLlJX1G36EaTRmQUhI7fITBh35IZNcuai7eTOMtt+BdeZ5ewCZANBrVi/84EULQ0tJCT09PyZ/RAkCjKQOxw0c4ceutpEMhAKJ79jL0yCN0PPggvlUrqzy6mYle/MfPeL8z7QTWaCaJTCYZfOiHpEMhnlpyMW+0GOp3OhRi8KGHkMlklUeo0RRGCwCNZpKkIxEiu3YB8K/rfoufL7/cei/yyi7S2pE566mtra32ECaEFgAazSRx1NRQc/FmACIuL2G3z3qvZtPFOHy+Yh/VaKqKFgAazSQRLheNt9yCIxAg6vQQdnkBcAQCNN5yi44GmoHcdddd/OM//qP1+xe+8AW++MUvcu2113LxxRdz4YUX8tOf/jTvc08//TQ33nij9fsnP/lJfvCDHwCwc+dO3vKWt7B582auv/56zp07V/F5jIW+MzWaMuBdeR6LHvhPkt8/SKS2kaYPfMCKAtJMji/+fA97zw6X9ZznL6zn8zetL/r+rbfeymc+8xk+8YlPAPDQQw/x+OOPc+edd1JfX09vby9bt27lne98Z0mO10QiwR/90R/x05/+lLa2Nn74wx9y9913c99995VtThNBCwCNpgwIIUh3LAMOEm1rZ/6f3ap3/jOYTZs20d3dzdmzZ+np6aGpqYkFCxZw55138uyzz+JwODhz5gxdXV20t7ePeb4DBw7wxhtv8La3vQ2AVCrFggULKj2NMdF3qEZTJqLxFAAj8ZRe/MvIaDv1SvK7v/u7PPzww3R2dnLrrbdy//3309PTw86dO3G73SxbtiwvU9nlcpFOp63f1ftSStavX88LL7wwpXMYC+0D0GjKRFgJgFgSKWWVR6OZLLfeeisPPvggDz/8ML/7u7/L0NAQ8+bNw+128+tf/5oTJ/IrLHd0dLB3715isRhDQ0M89dRTAKxZs4aenh5LACQSCfbs2TOl8ylE1bYpQggf8CzgNcfxsJTy89Uaj0YzWZQAkNL4OeDVWsBMZv369QSDQRYtWsSCBQt4//vfz0033cSWLVu46KKLWLt2bd5nlixZwi233MKGDRtYtWoVmzZtAsDj8fDwww/zqU99iqGhIZLJJJ/5zGdYv7462o2imndoDHirlHJECOEGnhdCPCalfLGKY9JoJkwkkUn4GokltQCYBbz++uvWz62trUVNOCMjI9bPX//61/n617+ed8xFF13Es88+W/5BToKqmYCkgfrW3OY/rTdrZixKAwAIRhNVHIlGUxpV9QEIIZxCiFeBbuAJKeX2AsfcIYTYIYTYMZ4iRxrNVJMtAHT5B830p6oCQEqZklJeBCwGLhFCXFDgmO9KKbdIKbe0teX1NNZopg0RmwAYiWkBoJn+TIsoICnlIPA0cEN1R6LRTBy7BjCiNQDNDKBqAkAI0SaEaDR/rgGuA/ZXazwazWQJxzOLflBrAJoZQDXDFBYA/yqEcGIIooeklNuqOB6NZlJEtAagmWFUTQBIKXcDm6p1fY2m3IQTKRwC0lI7gTUzAx2orNGUiYiZ/JVMSUZiOgxUM/2ZFk5gjWY2EI4n8Xuc1PpcOgpohvO5z32Ov//7v7d+v/vuu7nnnnuqOKLKoDUAjaZMhOMp/B4XAm0CKiuP/Tl0vj72ceOh/UL4ra8WffsjH/kI7373u/n0pz9NOp3mwQcf5KWXXirvGKYBWgBoNGUiEk9R43bidgqtAcxwli1bRktLC6+88gpdXV1s2rSJlpaWag+r7GgBoNGUCUMDcOJ1O7QGUE5G2alXko9+9KP84Ac/oLOzk9tvv70qY6g02geg0ZSJcCJFjcdJrdelw0BnAe9617t4/PHHefnll7n++uurPZyKoDUAjaZMROJJ2uu91Hrd2gQ0C/B4PFxzzTU0NjbidDqrPZyKoAWARlMmlBO4zufS1UBnAel0mhdffJEf/ehH1R5KxdAmII2mTETihgmozgwD1V3BZi579+5l5cqVXHvttaxataraw6kYWgPQaMpEOJ7C7zZ8AGndFWxGc/7553P06NFqD6PiaA1AoykD6bQkkkhZiWCgS0Jrpj9aAGg0ZSCWTANQ43FRa+76dSioZrqjBYBGUwZUKWi/6QMArQFopj9aAGg0ZUA1gzGcwG5Al4TWTH+0ANBoykAkYQgAv5kIBuiKoJpJc/XVV7Njx46KnV8LAI2mDCgNwC4AhrUGMKXIZJJUMIhM6u+9VLQA0GjKgPIB1LhdGR+AFgBTgpSS6KHDdH3ta5y87Ta6vvZ1oocOTzoP4+abb2bz5s2sX7+e7373uwDU1tZy9913s3HjRrZu3UpXVxcAJ06c4Nprr2XDhg1ce+21nDx5EoDbbruNj3/841xzzTWsWLGCZ555httvv51169Zx2223Wdf6+Mc/zpYtW1i/fj2f//zn88byve99jzvvvNP6/d577+WP//iPJzU/0AJAoykLEZsGEPBqJ/BUEjt8hBO33srAv/8H0T17Gfj3f+fErbcSO3xkUue977772LlzJzt27OCee+6hr6+PUCjE1q1bee2117jqqqu49957AfjkJz/JBz/4QXbv3s373/9+PvWpT1nnGRgY4Fe/+hXf+ta3uOmmm7jzzjvZs2cPr7/+Oq+++ioAX/nKV9ixYwe7d+/mmWeeYffu3VljufXWW/nZz35GImGYFb///e/z4Q9/eFLzAy0ANJqyYDcBuZ0OfG6HFgBTgEwmGXzoh6RDoazX06EQgw89NClz0D333GPt9E+dOsWhQ4fweDzceOONAGzevJnjx48D8MILL/C+970PgA984AM8//zz1nluuukmhBBceOGFzJ8/nwsvvBCHw8H69eutzz/00ENcfPHFbNq0iT179rB3796ssQQCAd761reybds29u/fTyKR4MILL5zw3BQ6TVGjKQMRWxQQQJ3PrfMApoB0JEJk166C70Ve2UU6GsVZWzvu8z799NM8+eSTvPDCC/j9fq6++mqi0ShutxshBABOp5NkEQGjjgHwer0AOBwO62f1ezKZ5NixY/zN3/wNL7/8Mk1NTdx2221Eo9G8c370ox/lr/7qr1i7dm1Zdv+gNQCNpixkfACmAPDqgnBTgaOmhpqLNxd8r2bTxTh8vgmdd2hoiKamJvx+P/v37+fFF18c9fg3v/nNPPjggwDcf//9XHHFFSVfa3h4mEAgQENDA11dXTz22GMFj7v00ks5deoUDzzwAO9973tLn8woaA1AoykDYSsM1HikdF/gqUG4XDTecgtDjzySZQZyBAI03nILwjWxJe6GG27gO9/5Dhs2bGDNmjVs3bp11OPvuecebr/9dr7xjW/Q1tbG97///ZKvtXHjRjZt2sT69etZsWIFl19+edFjb7nlFl599VWamppKPv9oiJlUsXDLli2ykjGxGs1E+eYvD/DtXx/m6F+9AyEE77v3ReLJNA9//M3VHtqMZN++faxbt66kY6WUxA4fYfChh4i8souaTRfTeMsteFeel2WKmQ3ceOON3HnnnVx77bVFjyn03Qkhdkopt+QeWzUNQAixBPg3oB1IA9+VUv59tcaj0UyGsNkPWC04tV4XJ0PhKo9qbiCEwLdqJfPv+jPS0SgOn2/CO//pyuDgIJdccgkbN24cdfEfL9X8lpLAn0gpdwkh6oCdQognpJR7x/qgRjPdUP2AFdoJPPUIl2tCDt+ZQGNjIwcPHiz7eavmBJZSnpNS7jJ/DgL7gEXVGo9GMxki8aQVAQRYTWE0E2cmmaenC+P9zqZFFJAQYhmwCdhe4L07hBA7hBA7enp6pnxsGk0pGM1gMgp1wO1gJJogndCRQBPB5/PR19enhcA4kFLS19eHbxyRT1U3lAkhaoEfA5+RUg7nvi+l/C7wXTCcwFM8PI2mJCIJox2kckjywvOk5EJOfO2bLPi9352VDslKsnjxYk6fPo3e9I0Pn8/H4sWLSz6+qgJACOHGWPzvl1I+Us2xaDSTQfkAVFkC0bYBLvodzv7oEeI/eZiOBx/Et2pltYc5Y3C73Sxfvrzaw5j1VM0EJIzt0PeAfVLKv63WODSacmBEATmssgT+pJHJGXL7ylKWQKOpBNX0AVwOfAB4qxDiVfPfO6o4Ho1mwkTiSWocWGUJ6uNGUtKQJ2C8b5Yl0GimE1UzAUkpnwe0UVQzKwjHU/hrPNRcvJnonr00xkYAGPTVAZMrS6DRVIppEQWk0cx0IvEUfq+bxltuwREI0BQNAjDoqZ10WQKNplJoAaDRlIFIwnACe1eeR8eDD9LxnpsRUhK95HI6HnwQ78rzqj1EjSYPLQA0mkkST6ZJpiV+j9MqS7Doz/+MpoCH6EVvwrdqpQ4B1UxLtADQaCZJphdAxsQjXC5aar30huLVGpZGMyZaAGg0kyScMMI77bWAAFprvfSOaAGgmb5oAaDRTBJ7O0g7rXVe+kZi1RiSRlMSWgBoNJPEMgG5czUAj9YANNMaLQA0mkmS0QCywzxba72MxJJEzW5hGs10QwsAjWaSWP2A83wAHgB6gpM3A8WTad7+rWd4/I1zkz6XRqPQAkCjGQen+sPcv/1E1muRYj6AWi8AvWXwAxzvC3Gwa4TnDvVO+lwajUILAI1mHNz73FHu/skbdA9n6voUdQJbAmDyfoAj3UZpiSM9I5M+l0aj0AJAoxkH24/2A7C/M2i9NhA2Fvg6nzvr2NY6QwCUIxLosCkADneHJn0ujUahBYBGUyIDoTgHuoyFf39npnfR3rPDzKvz0hzwZB3fYv5eDhOQ2vn3jsQYCusuY5ryoAWARlMiLx03dv9CwP5zGQ3gjbNDXLCoIe94n9tJnddVHhNQTwiPy3hcD2szkKZMaAGg0ZTI9qP9eF0OLlvRwj7TBBSOJzncPVJQAIBhBuopoAF85dG9fO6/3ijpuum05EjPCFesbAW0H0BTPrQA0Mw5hiKJCdnltx/r4+KlTWxY3Mjh7iCJVJp954KkJVywsL7gZ1prPfTmhIFKKXlk1xn+50hpET2dw1HC8RRXr2nD43JYDmGNZrJoAaCZc3zmwVe47fsvj+szQ5EEe88Nc+mKZtYtqCORkhztCbHn7BBAcQ2g1ktfTkG4Y70h+kJxBku05asd/6p5daxoDWgNQFM2dIcKzZxiOJrg+cO9JNOSoUiChhr32B8CdhzvR0q4dHmL5ezd3znMG2eGaA54WNBQuNtXS62HF45mawAvm76EgXCcdFricIxeKlpFAJ03L8B5bbWW0NFoJovWADQWr5wcsLJaZyvPHuwhkZJICa+eGix4zL5zw8SS2eUbth/rx+N0sGlpIyvaAridgn3ngrxxZpgLFjUUrfffWutlMJwgkUpbr718fACAtIRgdOzv+0jPCPU+F221Xs6bV8vJ/nDe+DSaiaAFgAaAaCLFLf/8At/85cFqDwWAoXCCXScHyn7eJ/d20VDjxiFg14n88x/uDvKOe57jkV1nsl7ffrSPi5Y04nM7cTsdnNdWy2unBjnYFSxq/4dMMlifLRLo5eP9OM1df3947AihI90hzptXixCC89oCpCUc7w2XNF+NZjS0ANAAEIolSaQkP331LEnbbvVAZ5ChSGXjzl861s+DL53Meu2LP9/Drd99MWsskyWRSvOr/d287fz5rGmvLyhgfrTjNFLCyf5w1ufeODvMlmVN1mvrFtSz/VgfybQsav+H/HIQ3cNRTvSF2bqiGcgkkY3G4Z4RzmurBWDlPON/7QfQlAMtADRAppxB70iM3xzpA+DMYISb/uF5/unpI1nHfuB72/nHpw+X5bqvnx7iQ/e9xJ8/8joHzSSr/lCcbbvPEU+m6SpDITXFjuMDDEeTXLduPhcvbeTVk4Ok0tJ6P5lK88grxs6/ezhz3e5gjFRasrTZb722tr0O9dELFhYXAG11ZkE4UwDsMLWOt62bDxjJZaMxFEnQE4xZC/+K1lqEyPgFNJrJoAWABsgIAID/MhfBb//qMPFUmlO23XA6Ldl+tJ/XitjPx8PpgTC3/+vLNAc81LidfOcZQ9A8tOMUcXPnf3YwMunrKJ7c14XH5eDKVa1s7mgiGEtyqDuT0PXsoR56gjFcDkF3MFPrp8us+zPf5uhdu8Aw+9T7XCxpril6zZZAtgnopWP9+NwOrljVBsDAGJFAR82dvtIAajxOFjXWaA1AUxaqKgCEEPcJIbqFEKVlxGgqhnL+Lmqs4fE3OjnYFeRHO04BcHYoswj3heLEU2kGQpMzC8WSKT7ygx1EEym+/+E38d5LlvKzV89a1TYXmoutXQCc7Avzj08fRkpZ7LRFkVLy5L4uLj+vhYDXxcVLDXPOTpsf4OGdp2kOeLhyVWuWBtA1ZAqAuowAWNdeBzCqAxgy9YB6LQ2gn01LmphXb7w+lgZgRQC1BazXzmur1RqApixUWwP4AXBDlcegIaMBvPeSJUQSKW7/wcs4hOCKla10DmV2w+rnUmzXo7H9aD8HuoL89bsvZPX8Oj565XIAPvYfOznVH+Ez160GDDOU4kc7T/H1xw8ULa0QiaeKRscc7BrhRF+Y6843TC8dLX5aAh52nRg05hOK8+Tebm6+aBGLmmoKawDmog3QVudl1bxarjR38sUIeJz43A56gzGGown2nh3mTcuaqPO6cDnEmN/j4Z4R3E6RZX5aOa+WIz0jpNPjF4QajZ2q5gFIKZ8VQiyr5hg0BkoAXLW6jf986RSnByJ86LIO6mvcvHC0j2QqjcvpsLSByQqA5w714HE6uHatsSAvbKzh5k2LeHjnaVprPdy8aRF//di+LA3geJ9hiuoPxWmr8+ad84P3bWdRYw1/d+umvPd++PIp3E7B9evbARBCcHFHE7tODiCl5N9fPEE8leY9WxbzxN4uBsIJ4sk0HpeDzuEYbqfIKvYmhOCJ9n+Cl3fCKDllAviNMw47IL5Dst0DDbvciNccfNF7CW+E7xr1e9p/Lsh5bbW4nJm92ur5tUQTaU4NhOloCYzy6TnIgcdh22cgXYYw2Y2/B2//y8mfZxoz7RPBhBB3AHcALF26tMqjmb0oE1DA6+KWLUv452eP8PGrV/LU/i5SaUnvSJz2Bp9NA0ggpRzV/DEazx3qZcuypqwuWh97ywoe2XWaW9+0FI/LwcLGGs4OZnbiJ/qMUsh9oRhQl3fOY70hXjs9xJejiazSzJF4iod3nuKGCxZYUTkAFy9t4om9XXz4By/z9IEerlrdxroF9ZZ/o2ckxqLGGrqHo8yr82XPNZ2GQ7+ABRthwUWjzvWNQz10D8dY3FRDR0sAT4MPjj7NFZHXeG4MU9qes8O8ZXW2lrF6vjH3A51BLQByOfEbCPXApg9M7jxHfw2Hn9ICoNpIKb8LfBdgy5YtWuetEPamJn94zXn8/taltNR6WdhgODjPDUVob/BZGkAqLRmOJkvOpLXTHYyyvzPIn92wJuv1lfPq+OWdV7G02VjUFjbWWA5oKSXHeg0B0F/Abi6lZDCcIJmWPLG3i3dfvNh6b9vuswxHk7z/0uwNxOYOww/w4tE+7n7HOj58+TIA5tcbtv6u4SiLGmvoHI7SnpvpGx0EmYYNvwdbPz7qfK+SRuJZVsbvtjtp3PXjrDwAKSVpiZUj0D0cpXckxvqcPAO7AHi7qdE8ubeLv3x0L//96SvzehPncnogzOIm/6jHzFjC/VA7H276u8md52efgoOPl2VI05lq+wA004RQzNAA/B4XLqeDFnOnrBa+c+bO3+4PGMuB+czBHn61vyvv9d8cNoqgXbky336+cl6dVfZ4UWON5QMYCCesrNlCAmAkliRp2sR//trZrPfu336S89oCXLq8Oev1LR1N/NW7LuSXn3kL/+uqFZaZRZmXlCO4aziaZf8HjIUGwN9SZPYZhBD55R78LdSlgwyFMt/nVx7dx3u+8z/W73vOGj0HcgVAwOtiabOf/V2ZCKZf7u3keF+Y3adHLxOx88QAV3zt1xyyfXZWEe4Df/PYx42Fv8U41wQCDmYSWgBogOJ9bRfkCIBzNpPMWH6Ab/xiP1977EDe688d6qXJ785b2HJZ2OgjGE0yHE1wvC/TCauvgBNYFVabX+/luUO9lnB648wQr54a5P2XduSZqxwOwfsuXcrSluzdsIrQ6TEdwV3DMUsrsAiblTwnutj4W3CQJhkatF569dQgu04Ocs7UsvaeMwTAugLf05r2Og7YupK9cnLQOsdonB4wNKozZQyvnVaEe0sSymPib4F0EqKzu+5StcNA/xN4AVgjhDgthPhINcczlwnFU3icDtzO7FuiocaNz+2g01yUzg5FWGYumKMJACklx3vDnOgPZYVtSil5/lAvl69sHbMI2sJG0/w0GLXs/6B8ANkoAfDeS5aSTEse39NJMpXmn545gs/t4HdsJqGxaAl4cQhj4R+JJRmJJQsIACNZbsKLjfk5R6TP+n5U9vFvDhvn3nN2iKXNfup9+Wa2NfPrONYbIpZMMRRJcMgMCx0rP0NldZdSg2gyfOnne/nFns6KXqMg4b7yCQB1vllMVQWAlPK9UsoFUkq3lHKxlPJ71RyPnc8+8rqVEDWTCceT/OEDu8ZMqIrEk1kOWYUQggUNNZwbipJOS7qGo5xv7khHywXoHYkzEksSTaTptmXzHuwaoTsY48pVrWOOXQmAM4NhjvWGcQhY0lxT0ASkhNHlK1tZ0RbgwZdO8r5/2c6ju8/x0StW0OAv3VfhdAja6rx0B6MFQ0ABmwAYex4FMTWHejlMMJYkmkhZ39Pzh3oAwwRUTEta015HKi053D1iLfrz671jCoDhKRAAUkru336Cn+WY4qaEcgmAgPl3Vaa+WYo2ARXhJ6+c5rlDpTXsAOOmf+Zgz7SLzd53bphHd59j+7HRdzKheIpAAQEA0F5vRP/0hmIkUpLzzSzY0TQAu8nmRF8mk/g5c3G7Yoz4eTB8AABnTA1gYWMNC+prCpqA1Fia/G5u2rCQ104P8caZIf72lo386fVr8o4fi3l1PrqDMZsAKLcGYCwwzSLIYCjB6QFDQPs9Tp4/3MdwNMGJvnBRAbC2PeMIfuXkIELA+y7p4OxQlO7haMHPgF0DqFx9p0giRSyZnUE+JaQShslmokLZjjLtaQ1g7hGOGzvXSKL0XdIrpwb50H0v8cLRid0w//bCcSvtv5z0mLvKsXZ8kXiqoAYAhh/g3FDUcgCvml+HyyEK7sQVKmIHsoXB9mP9LGvxW4v7aLTVenE7BWcHIxzvC7OsJUBzwFPwumpha/R7+P2tHfz+1qX8/I+uyIoGGg/z6rx0DY8hAFw14JlgNI0pOJpEkP5w3Fos37lxIb0jMX5qap/ri9QZWtYawON0cKAryCunBlg9r44rVhnnHM0PMBwx7oPhCgoAVd7i5FQLgIiZ1V0uJzBkfD2zFC0ACqB2mPb6OGOhHuCJJEjFk2n+4qd7KmJyKlUAhOJJAt7C4YPtDT66hqOcMXepCxtqaPR7Rq1jc6w3hMshcDoEJ20awN6zw6NWz7TjcAgj9HQwwom+kJG9W1tYAChzVEONm7Y6L39584VW/ZyJMK/eR08wSpcZCZQvAPonZ2owP9tMkIFwnFOmc/bWS4xQ1XufOwZgmdtycTsdrGgLsP+coQFc3NHI+oUNuByC104PFr3sVPgAlAN+MJyoqKDJY7JamR3tA5i7qAVmPAJA7RTDsfFnIKoInFiyfKWPFUoAjPUghuMpatxFNIDGGpJpyetnhszffTT53aOGgR7vDbGk2c/iphpLAxiOJjgzGGHdgtGjf+wsbKhh79lhBsMJlrcGaAl46A/Hs6p4giF467yuPCf2RJlXZ7RyPDMQodbrojZXOIZ6J7fT9PhJu3w0iyADIUMD8LocbFzcwIq2ACf7w7TWephXIONZsba9jheO9jEUSbBpSRM+t5O1C+pG1QCmQgDYBfSUmoFCKjKrDALAUwtOjxYAcxF1A0fGIQA6h4yFdiIdtULmZyoiAEZK0wDCo2gAC8zd766TA3icDpr9HpoCnlG1nWO9IZa1+OloCVg+gP3njLDF88chABY11lgRLh2mCUhKGMy59lAkQWNg/ElpxZhX70VKIxInzwEMxsIQmKStuabFEADhBCf7wyxp9iOE4MqVxnnPXzh6obk17fXEzXtm09JGADYubmT3qaGivii1EaikD8B+X0ypACinBiCE4UvQAmDu0WdpAKUv5koDCI1DaCjClgZQ/jZ/pZqAwqP4AFQy2O7TQ7Q3+HA4BM3+4gJASsmJvjDLWgN0NPs53meEgu4z49rXLsgv41CMhTZfwbIWP81mglquGWggHKexxkO5UJU/95wdzjf/QFmiTUSgxTABheKc6o+wpMmY6+WmABgrT0I5guu8LsvcddGSRoKxJEd7C/uTlAagfAGVwK4ZTqkfoJwCQJ1HRwHNPfrNOPNoovQdeacyAU1AA1CfiY3jeqVimYDG6OoVjhWPAlLJYOF4yhIGTQF3UR9A13CMSCLF8tYAHS1+gtEkg+EE+zuHafS7aS+0oBZBCQAhYEmzUcETMkJaMRBO0DiOUM+xUMlgsWS6iACYpA8AEIFW2pwjlhN4iVnx8/KVrVyyrJm3m5VLi7HaFAAXLW20ciouWtIIwKunCicwWQKgwk5gIaDO55piAVB6dnZJ+JszZqVZihYABZiIBqAiZEIT8AFYGkAZ2x8qMhrAWD6AZNEaMs0BDx7Ttq7q9Df5PQyE4gVr86sIoGUtAZaZxcqO94XYey7Iuvb6cRWQW9joM69bg8/ttCpy5moAQ+E4Tf7yaQDzbLX/8wRAKgGxockvNH7DBHSiL0QwlrRKPge8Lh762GVsWto06scXNvhYt6A+S1CsaKul1usqmA+QTktGzJIfFXUCh+M01LhZ3hrgZP8UZhyH+8BbD64y3QeqHMQsZtoXg6sG/eOMApJSWvXjp5MGIKUchw8glVcGQiGEEY1zsj9Mu1kcrsnvIZmWBGPJvExV5fRd3hogmjC+w2O9IQ50DvPeS8ZX0VWFi3aY2ccttVOjAbTWehDCKAVTvA7QJMMN/S00ErTq94y3QJsQgsc+fWXWa06HYN2COvZ3DucdH4wmkdLQpirpA+gPGcJ4SZPfMvtNCeWqA6SYAwJAawAFULvLWDKdF21S7PhEyjhuIj4ApTWU2wcwFElY4xpNAMSTaZJpWVQAQMYPoHbkTeZOfLBANvDx3hAep1HO2XBswrMHe4gm0uOKAAIjAgmwyh6rXX7fSCa72KhMmqCxjBqAy+mw2jnmmaysOkCTdAL7WwjIENGosXmwN32ZDGvb69l/LpinnSmzT3u9j2AsWbGkxcFwgia/myXNfk4PREp6hspCueoAKQKtRtXXVGXLZlQTLQAKYN9dRhJjL8qdtszLcGz8N4uKNoqXOQpImX/m1XlH3fEpDWS0MsLK9KMWw2Yz4kaVM/6Th17j7p+8Dhi7/aUtfpwOgc/tZEG9j6f2dwPjiwACqPW6+IOrVvA7Fy8CjPj3hhp3lgloOJJASiMLuJyoEMx55c4CVpi71UaM6KjReguPh7UL6gjGknkF35T9f3FTDVJmos/KjdIAljb7iafSVoBExSlXGQiFOldkYPTjZjBaABTAvriUYtJRN7jP7ZiYBlChMFAlAM5rqyUUTxXdiYWLVAK1o0w/yimrdtsD4TiJVJpHXz/L/dtP8sTeLo73hSzbP8BS0xHsdAhWzht/ctZn37GOLcsyqn1LwJMlpFU0UjlNQJBxBOf1AiibAMiUg2jyu7Oa2EyGte2GkFVht4qMADA0DbtWOJKzcZFS8j+Hey0T3ngYCMdpCngsgTZljuBwf3nKQCischCz1xGsBUAB+kNxq9FJKbkAKgdgWUtggj6AyiSCKfv/CrOh+EgRM5ClARTJAwBYNa8Wr8vBEnPxaFYCIBTnYFeQaCKN1+Xg7p+8zom+MMtbM+YMJQxWtAbwFUk2Gw/NAY/lpwEYtJWBKCcqFLSttlghuMk7gcEQAEvKZP4Bo1AckOcHGLZpAJARAEd7RtjwhV/wyK7T1rH/9eoZ3vcv29m2+9y4rz8QjtMc8FgmrSnLBaiED0Cdd5aiBUAOsWSKkVjSekhKcQR3DkcRQgmAiUQBGQ9ivMw+ANXQZIUZI14s9M/SAEZZnG/etIhn/vc1VlVNZYvvD8V5zQw5/NtbLqIvFCeWTLOsNaMBKPv9eO3/xcitB6SSwhon0J1sNG64oJ33X7rUalBjUUYnMEATQUuwloNas2HMvs5iGoBxb6v74WhPiLSEz/90D2cHI/SNxPjSz/cCcGSc9aki8RTRRJpGv5uFjTU4xBQJgHgYEuHKmIBmsQDQUUA5qIVlcVMNe84Ol+QD6BqK0lrrpaHGPaFSEJXUALwuB4tMx+2YAmAUE5DTrMujqPO5cDoEg+EEh7pGaPK7eceF7bxxdgX/9PSRrDo8KoJnPAlgo9FS62XXyYxdVtUBKmcYKMA1a+dxzdp5+W+EesHXAM5JChybBlBbRg0AjCSx/TkROPkmION3VYY6mkxx14930xLwEIwmafK7OW4r6lcKyifU7PfgNgMBpsQEVO4kMMiYk7QAmDuoQnDqISnJBGS2DPR7nRNyrIVjFRIAwRhtdV7LtlwsEqgUE1AuDoegye+mPxzntdODbFzSiBCCO69bzUVLGrnEZrO/cFEDHpeDy1aU5+FsCRiF6NJpicMhLBNQuQVAUcrWdMT4jpoJMr9MDmDF2gX1PLmvi2giZZndhqMJnA5h5TWo+0H5iu5+xzq+YO78P/XWlbxxdpjjfeNbvFUWsIoSW9Lkn8ECYPaXhNYmoBzsGgCUZgLqGo7SXu8j4HERjqeywu9CJUQFhSwTUKUEgLGwFxcAY2sAhWj0ezgzEOFgV5CNixsB8LgcXL++Pavb15JmP/u+dMOYiU2l0hzwkEpLa0c7GI7jMDNPp4RyCQCnG+mt54KmZMH+yJNhXXsdaQmHujImnKFIgnqfi/oa43saNu+H7mCUJr+bD715Gdetm8/5C+r5xDUr6Wjxc6IvVDDZrxiZvgyGAFja7J+aZLBKCACXFzx1ENICYM6gBIBKQCrk1B0KJ/jqY/ut9wwNwEeNx0kqLa2dfE8wxqYvPcHTB7pHvWakQrWAeoIx2mrtGkARE1BsYgKg2e9h+7E+0jJTgqAYzjHaP46H3GQwlXk6VovJslHGcEMRaOXty1x5fYkny1rT37LP5ggejiRpqHFbiXvKKdwTjDGvzocQgns/uJmfffJyfG4ny1sNn1aPraNbIe556hBvmJViVXkQFSa8tMVP70hsXIUVJ4Tyy0y2QF8u/matAcwl+kJjm4B+9toZvvPMEba9do5oIsVgOGFqAMYCqnbUZwYjxFNpdhwfPY7YHgY6nt3WWPSMlKoBjJ0HUIhGv9uql7RhcWk1/stBbjkII/Foisw/UJY6QBYVyjZd2uzH53ZkhYIORRLU17jxuhx4nA7rfug2NUUwsotdZtmPDquMR3ETTjSR4m+fOMj9208AGROQishSdaQ6K50LUAkNQJ1vLgsAIcTXSnltttAfimU5PAuZgF48auw2fr77bKZjVIPPsqErs48yURzoCuadw44SMlJiZe5OlkQqTX8oniMACmsAoQmagNRCvKS5hpbcUMkKkhEAxs50MJwYV8/fSVPOhKMKLTBOh2DN/OySEEORBA01boQQ1Plc1v1gaAD5f79lplYymiN40Nzx7zMFjRLKKiKryZYvUlHCfSAchnO+nARmd0noUjSAtxV47bfKPZDpgspiVIthbhSQlJIXj/bhdAh+c7iXPWeNB0z5ACAjNJQAODSGALAnj8XLVBBOObPb6rx4XU68LkdRDSAST+EQ4M0NdxwD5ehT9v+pQpVo6B3JmICmTAOIhyAZmfYCAIyEsH3nhi2tcjhqaACAKQCSRr0omwZgZ1FjDS6HyGrpmYta8A90BkmnJYOmOU5pEeoeGa15UFkI90JNEzgmn2eSxSwvCV30iRdCfFwI8TqwRgix2/bvGLB76oY4tfSNxGkJePC6HDhEvg/gcPcIfaE4H7psGWkJP/jNccDIFvV7jZtPmXSUADjRHx7VBmovHxGbQOZlIZTdViUx1fncRcNAQ/EkAY9rXFU6IVN6YSz7f7kpZAIqdxZwUcpec960MZfR9KdYu6COgXDCCvMcjiQs+7+6H4YjSeKpdEEB4HI6WNLstxr6FELlYEQSKU70h+k36wAp1M+jtQ8tC+UuA6Hwt8zZTOAHgJuAn5n/q3+bpZS/PwVjqwr9ISOLUQiB3+MiEs/ekb94zNgNfOjNHayeX8tLx43f59s1ANOpqpxsUhqCoxjhRMrafZcrFLRnxDBNqQe73ueyoj5yGa0h/Gg0mzvxDVOsAXhcDup9Ls4NGdElg2VuBjMqlWg6kowaSUxlRiXeKS1AmYAA6msMDUBVsc2rd2SyrMVvlfcuRL/NtLP/3DADobi164eMLyC3g1vZKadfxo6/2fjbxKe4wf0UUVQASCmHpJTHpZTvlVKesP0rmz4khLhBCHFACHFYCPHn5TpvKbx4tI//+1+v573eH4rTbEaZ1HicRBLJvM8taPCxtNnPTRsWGse5ndT7XJbZKJyjAQAcHMUMFI6lLBNGuUJBLQ2gTmkArqImoFA8VbQd5Gi8ff18Pnfj+WzuKE9453i4cnUb2147R38oTiieKnshuKIoAVCuaJMKJhuta1cCwCjVkUhJSwDUed0Eo4k8TTEXo6Vn8VBQu2ln37lhowyEzRxXbyYMTokPoCICwPz7RGanGahqUUBCCCfw/zD8CecD7xVCnD9V13/s9XP8x4sn84pg9YXiVtcpv8eZ5QSWUrL9aB9bV7QghODGjYYAmF/vRQhhLaKWDyCcoLXWaKZSTAAkUmniqbRlwiibBmA+2K02E1AxJ3AknizaEH406n1uPnLF8rKGeJbKHVeuIBhL8s/PHgGgMTBVGkC5u05VrtxAg9/NosYa9p4btjYjKgdAbQiUeWheob7HGBpAKJ6y6krlokw7HS1G6YmBUDyrJpMQgsaa4t3jykYlTUDq/LMQUc6ww3FdWIjLgC9IKa83f/8sgJTyr4t9ZsuWLXLHjh3jv9iL/wT7H816aX9nkP5QnE1LG63FLw28eKSPJc1+ljTV8OrpIXwuh9V7NZxI8erJQc6bV8t8c2f9+pkhXE4H69rrrJDPFW0B2ut97O8MEkmkEELgdTpYV6AUQjIteelYPw01boYiCTYsbqB2ArvxXI72hugdiVkZuQe6goTjKTYVsNfvOTuMBC4YowftdOONs8OMmHXtV8+vo7V2CoRA8Bz0HYa7jhtOx8lycjvc93ZovxB8jZM/Xw77OoNEEynWzK/j1VODrG6vozXg4VhfmO7hKEuajJ7NlyxvxlVAkA+EE+w7N8wFixqoL5Bod6w3RHcwRqPfw0jM6D/RXu+zIogAXjk1iN/jZM388pQCKciJ38Dln4HrPl/m874A378B2jeUP8JovFz3RVi8eUIfFULslFJuyX29mqUgFgGnbL+fBi7NPUgIcQdwB8DSpePrJmUh05DOdq6mUgmcIk0ikaDG3PwmU2mcIo3HYRzvFilIZz4bDEdxijQNXmG9tm6+eaOnUziROEUamUpCOoVMJ/E6wOMSjETjeWMwPmZc0+swP5tOQnryO+pUMoHPKa1rukUa0smCY0AmcTkchd+bxiyqd3OwO4ZTqPlNwfgD82DJ1vIt1vPXw+rfguhQRcZf54FgOEEiYdzvbtS9nQaZIpFM4HakcZIqeN/VuIz7MhaPgyf//XQqic8pqfPAYCiJAOv5UXgdknSqyL1XLjouhzWlByd2DkX55i8P8OWbLxi9Qm37BRX9+4yPCmzWpZRV+Qe8B/gX2+8fAP5htM9s3rxZlotLvvKE7Lhrm/yvV05br+07NyQ77tomt712Vkop5fvufUG+6/89b73/ift3yq1/9aRMp9MFz5lMpWXHXdvkt544IKWU8h1//6z88Pdfkt/+1SHZcdc2ORyJ533mcHdQdty1Tf6fR3bLjru2yecP9Ux6bqlUWr7pL5+QH/v3HdZrX/r5Hrnuc48VPP66bz6ddexMIZVKy2u+8WvZcdc2+frpwWoPZ1ry37vPyo67tslv/vKA7Lhrm3z15ICUUsp/ee6o7Lhrm7ztvu3yiq89VfTz8WRKrvjso/Lrj+8r+P4HvrddvvMfnpNP7euUHXdtkx13bZMPbD+RdcxH//Vlef23ninbnMrBT3adlh13bZOvnRqo9lCmBGCHLLCmVjMT+DSwxPb7YuDsVFw4mUpbNnJ7tyJVY16FGda4XVk+gH1nh7nILHpWCKMDliMrD6Chxs1qU/U9VCASSEUMKSdwOcpBvHZ6kO5gjLevzzQLr/MZc0kUyDMITzAKqNo4HII/vGYlHpfDalSjyeZ806z34lHDht1gywMAONITYl5d4QggMDqwLWmqKZoNPGg2f7GX+s7NyWjyuyvvBB4nKtii4r6JaU41BcDLwCohxHIhhAe4FSPktOL0jsRRzbG6hjPOLVUGQtWa8XucWYlgfaF4wYxJOwGPKysTuKHGbdk+D3bmO4JVxJAKnZtIY3gpZVZ/11/u7cLlELx1TUYAqPjvQk1hwmYewEzkdzYv5pXPvc0S2ppsljT5CXicvHpyEMBKBFP2/FMD4aIRQIqOlkDRbGCVONle77OES25EVpPfqN4qS/Q3/t2TB/niz/eUdOxEiZkboYqHp05zqiYApJRJ4JPAL4B9wENSysr+1U1U/DhkawCq0XizLQpIJXAlUmmGIgkr9r0Yfq8ROZRKS4LRJPU1bhY31VDjdnKwq4AGEFcagPHQTCQT+JMPvMLH/mOn9YD9ck8nW1e0ZJVHGK0eUDieGncZiOnEREJY5woOh2DtgnrrvlILv9oQSFk8AkjR0eLnZBENYMAUAEIIK1giVxg3+j3Ek+mSemtIKXnwpVM8f6iyyVdKA+gvIUP5uUM9FW9sf6o/TG+RSKtKUtVicFLK/5ZSrpZSniel/MpUXVct+nU+l9U1C+D0QASf22HFMdfYBIBSYVWVw2IoDUCFXKoqlavm1xYMBc0IgIlrAK+cHOCXe7v41f5uDnePcKQnlGX+MeZqVoDMCQVV1UvHWwhOM3NQ0We1XpdVosHef3gsrXZJk59gLJmV1wKGuTIUT1nPhDIDNQXyTUBQmrnlZH+YzuFoSWXUJ0OpJqC9Z4f5wPde4pmDo1f0nSyffGAXX962t6LXKMScrAbaOWQIgA2LG+gKZjSA431hljb7rbLCfo+TcMKo76+6To2pAZhmI/WwKLV4ZVttwfZ6qmxEJg9gfD6AVFrSZfoz/vLRfTxq9nB92/nZAqC+iAaQqQQ6czUAzeicv8AIX7SHcdp7JxQqA2FH9cY4PZCtBahCcCru/12bFvGezYuzEsFgfPWAtpuZ9qEKl49WAmAsE1DnsGEt6A+Vz1fQHcwXcP3hOOcGK1wxtQBzUgCcG47idgrWttfTPRyzTCcn+0NWCVwwMnxTaUk8laYvlG0eKkbAa2gAuQKgtc5LXyieZweN5GoA40wE6wnGSKUl162bz7HeEN/+9SE2Lm5gQUO2U7RYTwCrGYxXC4DZitIA6mvyTYLAqE5gyJRGPz2Q3dgloxWbRQGXNPKN92zM68swnoqg281Ku4X6cJQTtdEaSwNQBQfLqZHc+s8vcs9Th7JeiyYya8xUMicFQNeQ0cClvd5HJJEiaCYTnewP02HrzVpjmkUi8ZRlKxxLAKjs4VwB0Oh3F7SDhnKdwDYB8MD2k7x2anDU6501/Rnvu3QJV69pI5GSvH19e95xxXwAE+0Gppk5rGmvQ4jMvQjZJqDSNYBsAWCVfh6jDMd4TEAvHTeilRIpWbAsipSSR3efY2iS0TuWCWgMrUTNMbdiwGQ4MxjJa7ITTaSqEpE0JwVAp9nCUTm/uoejdAdjRBNpOlozGoDf1uBloEQBEPC4CMXzNYDMLihnBx4zSjHXmY5MuwD4yqN7+c4zR0a9nlIbFzTU8Pmb1nPJsmZ++6KFeccV6wmgdjY1bu0DmK34PS7WzK+zelyAUVDP5zYe/7F8AI1+NwGPM88ElDGLjv5MlFoQ7uxghFP9EZaY/ZEL7brv336SP3xgFz/fPbmIceUUH0srUYEh5RIA8WSaWDJNNMfUG0umGQzHK+5szmVOCoCu4RjzG3xWc+yu4RgnzJrndg3ALgBUiOhYRcf8XifhWL4GoD6X+xCE4ykCHhcOh8Djclg7k3RaEoqneMUM3yuGimha0OBjeWuAhz52maWy28k4gbNvZKWRBLQJaFZz321v4i9uzC61VedzI8TYC7gQgsVN/uImoDF6MSgNYWAMO/pLpv3/2rWG/yqUYwY60jPCXz5qOEqLFTYslYwPYPQxqb4ahcKnJ4LagEVtwR7ptKHtpGWmgvBUMecEgJSSc0MRFtTbBUDUqnm+LMcHABkTkL3RRTGKaQCZXVCuDT5pJWF5nQ7LNqkW5s7haFbY6munBjnVn9mJnRuKUuN2Zqn3hfC4HGZTmMIagDYBzW4WNuZ3bav3uWgJeMe8p8EwA+UJgJz2j8VwOx3UeV1j7ra3H+unzufiYrO6rD0JM5FKc+cPX8XndiIK9OkYL6WGgfaGyusDUIIrajMF20O/+6c4L2HOCYDhSJJoIk17g89SfbuGY5zoD+FyCBY2ZtRkFRoZSRgCoKWEZCO/x0U0kWYgFMfjzKjZxRxhYVspZq/bYZmA7DfcrhODgPEQfPC+l6xdEBgawIJGX0nNXIyKoDkagOUD0CaguUadzz2m/V9hCIAcE1A4Qa3XhaeETnJNAU8JAqCPNy1rtsyV9mfge88fY/fpIb767gupcWdX6Z0IKhEskkhlLca5qLajwQoKAPvPFe+clsOcEwCqOfX8eh8Br4s6r4uu4SjH+8IsbqrJ2g3V2Or7q0YxY6F20ueGotSb/VehuCMsHE9an/G6nFYegN3m+MpJo6n8juMDDEUSHLBlFJ8djLKwobQyCKoJCBihsFLKCfcD1sx8br5oIe/ZvLikYxc3+QlGs3MBBsLxkjuxGeUgips3eoIxjvaEuHR5s5WVHoplFsbXTw+xojXADRcswO9xjUsADIUT/OSV01mv2fNtRjMD9ZU5CqiQCcj+cymJaeVkzgkAZU5RDrG2ei89wRgn+8IstZl/IN8ElJvgUggVTnluKEpDTWZXrbJyB3P+wKFYJgvX43JY6qD9Bn/FjAR6al8XkN1isnMoyoKG0cP4FKoN4L+9cJytf/0UX9q2l4iVB6A1gLnGbZcv5/Yrlpd0bKFcgIFwaZsiMMxEozmBVfP6DYsbrefB7gOw9zM2Iu1KX5D/8+WT3PnD17Iib+xml2KaiZQy4wMokwBQPjh7vo/9Zy0AKozKAm437f/z63x0Dkc53hfKqmEO+U7gUkxAavfSORTNsst7XU78Hme+BpBIWYuv1+WwegKrG+78BfW8fmaIeDLNk/u68LkdSGk4xJKpNN3B0gVAvc/Fy8f7+Yuf7mF+vZfv/+Y4P3vtbNZcNZpCFMoFUGUgSmGsgnCqDML8eq/VD8O+yAejScs0lNuoaSwOmSVYsuzuyRTKalrM7BKMJS1BUS4BMKYGoH0A5SeeTFtmk84hdaOZAqDey6GuIMFokqXNxQRAkoFxmoC6hqN5jtkmv4fBSI4PIGY3AeX7AK5c1Uo8mWbb7rMc7wtz65uMnggHOoN0BWOkJSwosRJmnc/wT9x80UKe/tNr2Li4gZePG+aliXQE08wdCuUC9IfjJbfibPR7GBwlCsjqYFfntbToEZsJKBjNNLQfrwZwuKeQAEhb3fKKmaZUdWCvy1HGKKB8DUD7ACrMn/94N+//lxcZCifoHI7SEvBYjqv59T5LLVuWawIyF+buYIxkWpYkAJRDN5mWeQKg0e8uEAWU0QDsYaDKNn/FKqMn6beePAjAR65YbrSY7A5ybjATAloKv/empfzp21fzt7dcRI3HybffdzH1Phc1bmde9qZGY6dQLsBgKFGSWRSMzU8wliza87p3JI7XZUQLKS06HMvWAJRmEPBm+wCklPztLw+w79xw3nmllBwxy7Dbc2ziqTTz65UAKLzoqszcpc3+CkQBZcZiH1c5S06UwpwQALdfsZz+UJy/+eUBOociWQkx8+ozP3fkmYCMG07tesajAUB26j2osri5UUBJKwbf63JaOwN1w62cV8v8ei+n+iOcv6CeJc1+VrQFONgZ5JxZ06jUWvhvWd3GJ9+6ylrslzT7ufeDW/jT69eU9HnN3CU3FyCeTBOMJUs2AamCcbkasKI3GKO11uitrbRRez0guwmoxu20+miAEclzz68O89COU+TSNRyzzDe5GsB8swRGsV23KgPRYfZFTpchSStjAiqiAWgTUPm5YFEDH7xsGf+x/QSvnBq07P+AtQsQwlgQ7TjN5Cy16xmPBgCUrAHUjGICCnhdbFpixEVfZxZ4Wz2/joNdI3kO7Ylw6YoWPlKiI1Azt7HnAqiFvFQNoFgejKJnJEarGZLqcAgCHqelASRTRgkVlcwY8LoIJzI7cmWeOVyg4ZL9NftOO5ZME/C6CBTwyymUA3hps2EZyE1MmwhKA0imJUnTv6AEQJPfrZ3AleJP3r6atlovg+EE8xvsAsD4ub3eV7A3qN/j5Ix507eMUQlUHa8oJADsEj6ZMtLClcrrddtMQOYOJ+BxcXFHIwDXrZsHwOr5tZwZjHCoa4Rar8uyjWo0lcSeC6Cyekv1AVh5MEUWuJ5gjLbajDDxe13Wgqt28JYG4MnWAFSMfmEBkAmZttvd48k0HpeDpkDx6KR+ywRUkzWOyRCMZYRN1HzWlWBa0FCjNYBKUedz83/NVPgFNg1AJYPlmn8UfrfTyh1oGqMXAJDVWauQCWgokrBUyXAiOwbf47RpAPEkPrcDp0Pw/ks7+M7vb2bD4kYAq8Xkc4d6S7b/azSTxZ4LYBVHLNEE1DhGQbjekXhWUlrA47Q2QWrXrARAICcKSGnL54aieYv0YVsJdrvdPZ5M43U5Cppl7WOq87loNp3F5fAD2BMx1c5f/b+w0TflGsCcCv6+acMCYokUV61us15TpXA7mgMFP+PzOK32kSVpAN7RNAAPUhoxzY1+T14Wrt0HMBLLdnrdcEGmwqcSAJ3DUVa3t6HRTAX2XAC1ay7ZCRwoXhAulZb0h2JWVA5gJnsZi6VqYpTRAFxEEoZN3uEQWRE6R7pH2Lik0fr9cPcIzQEP/aF4QQ2g0e+mv5gJKBSntdZrFWqcbP0hYy6Zc8RyNID2Bh/BaNIa21QwZzQAMBxZ79myxDL7gKFOfvCyDt5ZoIImZHbnNW5nSY3TPU5j1w75AiA3Gzi3Do/dBGSEhxaWz0ua/VaJiYVaA9BMESoX4GRf2IpXH08eABSOcx8IGz267QKg1usqoAGYPgDzeYnk5MwAHMoxAx3uDrF+odGpzJ79G0uZJqBREtT6RmI0BzyWX8+emTxRgtGElX+QqwGoHh5T2ad4TgmAYnzpty/g8pWtBd/zm2WSS814FEJYC3qhPADIePpza/HbTUAjsVTRXrdOh2DlvFpgcg5gjWY8LG324xDwR//5Cn//pNHQpNRSEDVuJ16Xo6AT2MoBsGsAXmfGB5BjArInaEK2ALD7AYbCCXpHYqxfaHREUwutlEb1Ta/TQXPAU9Qvoep/KU18JDb5EM1gNGmZzdR4Mj4A41meymQwLQDGQO36SxUAkPEDFHICQ0bCqxu4WDG42lFKNK+eZ5iBSq0DpNFMlga/m5984nLuuGoFbXVeNnc0FQycKIQQwrC3F1hsVRZwtg/AZWnIynFaZyWCZWcKq+OaA54sAXC4x3AAn680APPZUtm9ygQ0HE1aETnZ44rTUuu1CYDyaABqnsonEUsYWclqMzeVfoA55QOYCP4JCADlByiqAZgRFGqHkwkDNVpQJlNpwvHkqGV2V7cbAmBBo9YANFPHxiWNbFzSyJ/dsHbcn13a4udAVzDvdSUAWu1RQDZHb64TOFcDUFFAFy1pzOq7rYSBMgGpBVeZWT0uB3Uu41yDkUSWBpI2/RItAQ+15nVHopPTABKpNNFEmrY6L/s7g1bZl6jpkFY+xrH6JpQTrQGMwUQ1AJdD5NXXsTQAs6KicgIHbJnAYOxQ7E7gQlyxspUFDT7WtteXPC6NpppsXdHCG2eGLKeuojdo7Hhb7RqA16YBmAJAPQ/2Kr1gmIhcDsH6hfWc6AtZppXD3SN4XA6WtQRwO4XlBLYEgNNR1Dk9GEmQltBS67ESNSfbqF7No80UNKorWDSRwutyWlGGuSagcDzJ7tODeb08yoEWAGMwIQ3AYzRoya3RX+9z4xCZmy3PCWwKgFgibfYJKK5eX7CogRc+e23J9dw1mmqzdUUzaQkvm52/FD0jMWM3btvwBLxOQvEUUkqGowmzt4bTfE+ZgDJZ8wGvi5XzaklLOG529zvcPcKK1gBOhzAj7LJNQF63M+Ocztl1qxyAllovXpcTj9Mx6SggtYDnm4DS+NyOorkS+84Feee3f8MOs25XOamKABBCvEcIsUcIkRZCbKnGGEpF2RvHpQF4XQU7dDkcgoaaTDJYJCcPwGuqo7GkoQHoEs2a2cTFS5vwOB28eLQv6/XeYIw2swyEwu9xkUpLYsl0VhkIyBQutKKETG1ZBUYc7h5BSsnBrhHOM1/zuR2WZpClARRp1KTKQLSaz33A68zKAzjVHx53aQglQJSpyYoCSqbwuZ1G5zSfK88HoNpE5uYVlYNqaQBvAO8Gnq3S9UtG7TrGIwA+eFkHn3zryoLvGYknKgw0xwmsNIBkynQCawGgmT343E42LW3kxaP5GkBrjiYbsNn5R3IEgHpeIomME7jW6+K8tlqEMATAD18+xZnBCFea0X1ZGoDNB1DMBKTKQDSbfolan8uKNuoejnLN3zw97sb0wzkagJUHkEhbz77KWbCT2162nFRFAEgp90kpD1Tj2uNlIiagq9fM490XF+60ZNQDMjWAeBIhMgu/8gEMR5KkJUXDQDWamcrWFS3sOTuU1VmsdySeVQYCsMXeJwlGE1YEEGSeSbWBGoklqfW58LmdLG6q4blDvXx5214uW9HCLVuWANl1tmJ2AVAkQ9kyAZmO2YAnIwBO9IdJpiWvmo2acnloxyl+6++fy9MQLB9AXWENAAoXjJx1AmA8CCHuEELsEELs6OnpmfLrq5utlGYwpWCEwhl/0CM9IdrrM/18lSBQTqDRfAAazUxk64qWPD9ATzA7CxhsAiCezDMBqWdSBVGM2EpFr5pXx84TAwgh+MZ7NliVb71uZ17cvcfloMbtxONy5Nnde0fiCJFJYKvzuax8hLNmGfaDBSKaAH6y6wz7zg1zZjCS9Xq+AMgUg5t1GoAQ4kkhxBsF/v32eM4jpfyulHKLlHJLW9vUlz1QX7o9e3gyNJr1gBKpNM8e6uEttrIUXnMXoHYfAe0D0MwyNi1txOPK+AEKlYGA7F2+vReA8V5GOEB22ZRVps3/L24638pchmwNQJmAvE6HmZ+Q362sLxSjscZt9QgP2IrTdZpl2A925Refi8RT7DxhOGvtvbsh4wTO9QHEkulsDaCAAPB7nBUpD1GxFUZKeV2lzj2V3HBBO//+kUvySkVPFFURdNeJAYLRJFevmWe95zFvNhWRoDUAzWzD53Zy8dJGXjxmCABVBiI3mi0T6ZNvAjKiehwZDcAmAH5/awdLW/x5ze7t7VYzUUDG82b3yyn6zCQwRa3Xxck+oxKq6sPRE4zRn9MpcPuxPuv8B7qCVgl3yM5n8LgctjDQNC0B09JQ68kLAx2KJCqy+4cZYAKqNl6XkytXlU/zaPK7CcdT/GJPFy6H4PKVLZlruZUAMDUA7QPQzEIMP8CwVaoByDcBqV2+qQHYTUCQnSgWspVNWdLs5/2XduSFYPvcTqv8ciYKyFh0633uvBj7gXA8q9JprTfjA1B9OCDfDPT8oV48Lgdtdd6CGkCNGe3jczms2kSxRCpLGEUTaUu4wSwUAEKIdwkhTgOXAY8KIX5RjXFUA5Xd+7PXzvKmZc1ZOxvLB2BpAFoAaGYfb1ndhpSw7fWzmSSwPCew6gucZCSepD5PABgmmXRaWk7g0cjSAGw+AMiO8FGMxLKFjl0AdA5FWWNW5M0TAId7edOyJi5YWJ/3nl2Q+XJ8Ej6XCjbJTwYbiiQqEgIK1YsC+omUcrGU0iulnC+lvL4a46gGKu64dyTGNWuzNQslAJQNUPsANLORi5Y0sn5hPf/2PyfoGTHMKblhoMrO3x2MIiVZGyXjfSeReMrqqTFa3SwwFtxMIpjxGUsAeF15Td9HotlCpdZn9CFOpSVnh6JctKSRep8ra5ffHYyyvzPI5StbWdNez5GeERK2GkPFBEA0kbKq+6r1QTWkByMPYFZpAHMZewela2z2f8gkgqkoAO0D0MxGhBB86LJlHOgK8ujuc0AhH4Bx73eZ9vbcHb7RMSxlLdy13tEXSLsGoEwvY2kAdsez+nkoYpitFjT6WNNel7XL/83hXgCuXNnGmvZaEinJ8d6Q9f6wzZdhJKZlwlLVs6/8Dr2mGVhdUwuAWYIyAS1qrLEyFxW5YaA6EUwzW3nnRQtp9Lt5cl93XhkIMLJ9hcDqxpfnA3A7icSTVonmMU1Atkq7VjVQM+iizpsvAIK5GoA5vqM9I0hplG5eNb+OA51BpDTi/Z871EuT3836hfVW0yZ78Tu7BmBv/mTXAJQprDeoBcCsRBV8umZtW56jypNjAtKlIDSzFZ/bye+ZSVq5ZSDA0BICHhedw8ZCmGsCMkozpKwSzWOagFzO/FIQNhNQNJG2zDXxpNGrO7s2kfGzajjT3lDDmvl1DEeTdAdjpNOS3xzu5c0rW3E4BOe11eIQcNBmIhqOJqz+3UoDSKbSJNPS0gCUM7zPXAPiSaMumBYAs4T5dT7uuGoFH758ed576iYYCMfxOB1T1hZOo6kGv7+1AyHyHcAKv8dJpxlxk6sBqLaQJZuA3PmZwF6bCQgyxRmVNpBlAjKPOWTG/i9s8GV2+Z1BHt55mq7hGDesN1q3+txOlrUGimoARlRSyhqL0gD8Hic+t4M+MzpKJYGV2nhnvOgt5hTjcAj+zzvWFXxPLfhGGQht/9fMbpY0+/nwm5dn+cXsBLwuTpiVPfOigNxGcTa1WI/1vHhdTpJmrw17MTjILPTBqNGDwxIqNq2j1tIAjAW9vcFnxf9vP9bHA9tPcsmyZm7csMD6zJr5dezvtAuARJYJqHckbmklKhFMCEFrrdcqRlfJLGDQAmBa4XQI3E5BIiW1+UczJ/iLm84v+l7A60SV08nd4fu9RhSQEgB1Y2gAaocdS6aJp9K4ncIqE6EW5ZGcDmSFnMCHu0eo9bosk1RrrZfvPHMUgC/ffEGWKWv1/Doe39NJNJHC6RBEE+ksJ3AskbJyE7w2bb+l1mvlRwxVsBIoaBPQtCN3V6LRzFXsm6BCiWCheNLq0lWKBgCGAIgl0tZzBhnhogSA0gDqCziBzw1Frd69AGvaa0mlJR+5YjlrzC59mffqkNIQGrldzVQYaCxHAwBoq/VYGsBwhTUALQCmGV6r6YU2AWnmNqoktLNAdz2/x0VaZkKmx4oCUhpANJEinkpl+dcyLR9zfAAFBABkevcCbF3ewrIWP5++dlXeNe0+ApVpnBUGmkxboaBqfGBUIM31AWgT0BxBqYI6C1gz1/Gbz0Ct15UXJaQEQncwhsfpsHb4xbBrAPFkOlsAKB/AKE5g+/No1wD+6NpVfOKalTgd2eMDWNbix+9x8uS+Lks7sDQAMypJhYLax99a56EvFCedlhUXAFoDmGaoG1NnAWvmOrXmM5Br/oHM89EdjJWkLdubLeUKgLocDcDqQWy7rseVicpb0FCTde5Ciz+Ay+ngD646j8fe6OQXezqzrqWikpQG4M3RAFLm4q8FwBxD3ah+bQLSzHHUM1DIH1ZjaQDRMc0/kLGxRxOGE9i+41bnV0llxRzL6ji7BjAWf/CWFSxpruGfnj4CkMkDcDlJpaUVemr3AaiyGL0jMasUtNtZmaVaC4BphroxtRNYM9dRu/x6X/7uV+36u4djJWnLlgaQMDUA24Lq9xhZxxkNIIHTIbLs8mATAI3ZGsBo+NxOPn/jepJmOJPdCQwwaO7w7VFAqg9x70icoUiCxgrt/kELgGmHR/sANBogowEUMgHVuI3XekdiBd/PxWsLA43lmICEENR6XRkfgNmAJtfvMBENAODadfO4Zo1R+FGZcpRwUSae0TSASoWAgnYCTzssJ7BHm4A0cxu14Bb0AZjCodTe2UqzNhyv6bws+zpbRdBgTiG43PG0j1MAGO0pN/I/R/qsWmAq2m/IrPuVlQdgagB9pgColP0ftAYw7dBRQBqNgd9yAucvgPaw0FLMpVmJYMl01oILxvNmzwMoJHRqfS5qva6CJqmxaK318s6NC23jyTYB2TWAJr8HhzBMQJUsBQ1aA5h2qJ2KFgCauY7SggstxvYksVIEQF4YaI5T1V4SOrcZjKKjxZ/XN2Ci+FzFTUAOh6A54KUvVHkNQK8y0wwdBqrRGKhNUKEon/FqAN6sRLB8E1Ct12WFf47EkpYZxs7d71hHyiz9PFksDSCc7wQGo0BeTzCuTUBzjYwJSPsANHObgOUELmQCyq/UORq5GkDuglvnyzYB1Ra4pquEhLNS8do0AKMGWK4A8NI5HKloKWjQAmDaoXYqOgxUM9dRtvZCYZAelwOXmYA1Hh+Ayr4tpAGM5QQuJ0oDGI4kLHOQndZaD8d6jEqoDRUqBQ1aAEw7PE7jxtDVQDVznZXzavm737uIt50/v+D7KhmslMVa2fwLlYIwzuEe0wlcTuxOYK87X6toqfUSihtlIrQGMIfQGoBGYyCE4OZNi7IcpHaUn6yUgAkhhNEXWJWCcGafUzmB48k0kURqCjSAjAmosAaQ6ZGsBcAcQpeC0GhKQzmCS/EBgLHrjpmlIArlAQB0mT2Ip8oElErLIhpAxgmtBcAcYn69jzqvq+IqqEYz0xmtVlAhvC4HkXiKRErmm4DM561YE/py47M5k3Md0mD0SVbMujBQIcQ3gJuAOHAE+LCUcrAaY5luvGfzYt52/vyyRRtoNLMVvztTLroUfG6nZefPXXTtDV+g8gLAXv2zkIlrtmsATwAXSCk3AAeBz1ZpHNMOl9ORZf/TaDSFmYgGMBwtHHdvaQBmE/qxmsxPFq/LgSo1VEgDsK8BlawFVBUBIKX8pZRSpdS9CCyuxjg0Gs3MxT+OKCAwdt2qxWIxH4DSAEr1K0wU5ZSGwhqAajgfqGApaJgePoDbgceKvSmEuEMIsUMIsaOnp2cKh6XRaKYz/nFEAYFhd1fZvoVKQQB0Dk2NExgyyWm5ZaeN15zU+VwVNf9ABQWAEOJJIcQbBf79tu2Yu4EkcH+x80gpvyul3CKl3NLW1lap4Wo0mhlGwOPM6tQ1Fl53xgRUKBEMps4HAJmFv5i/r7XWW1HzD1TQCSylvG6094UQHwJuBK6VskwFNjQazZzhPVuWsNJsvF4KXpeTYaUB5JmAjIV2KjUAZfoppAEALGn2W9nOlaJaUUA3AHcBb5FShqsxBo1GM7O5YFEDFyxqKPl4n9tBPGn04M01AVkdxoJRHCK72FylUKGgxTSAb75nI6Ky63/VqoF+G/ACT5hdd16UUn6sSmPRaDRzAPtCm5t85XI6qHE7iSRS1Pvyu4FVArXzL6YBtNVVPhqwKgJASrmyGtfVaDRzF/tCm6sBgOEIjiRSBauPVgKvZQKqXs7PdIgC0mg0mopj1wAKOY5VKOhU1eFSC3+hPICpQgsAjUYzJ7AvtIUWXRUKWukcAIVvlDyAqUILAI1GMyew2/0LaQC1U6wBqPEUKgY3VWgBoNFo5gT2XX9BH8AoLSgrgdIAtAlIo9FoKkyWABhFA6ibYh+ANgFpNBpNhbEvtKP6AKZMAGgNQKPRaKaEkjWAKQoD1RqARqPRTBG+sZzAUx0FpASA1gA0Go2msozlBK6bYh+AGo+OAtJoNJoKoxZaj9NRsNRD1TSAIqUgpgItADQazZxAmVqKlY9WXcCmygmsSj1Plc+hEFoAaDSaOYGlARQRAOsW1LG2vY417aWXmJ4M16+fz/0fvZRFjTVTcr1CVKsaqEaj0Uwp3jESrxY3+Xn8M1dN4XicXL6ydcquVwitAWg0mjmBbwwNYC6ivwmNRjMnUDv/QhFAcxX9TWg0mjmB1gDy0d+ERqOZE3jHiAKai+hvQqPRzAm0CSgf/U1oNJo5gcvpwOUQVc28nW5oAaDRaOYMXpdDawA29Deh0WjmDF63s6rll6cb+pvQaDRzBr/HWdXyy9MNnQms0WjmDF999wbaG7zVHsa0oSoCQAjxZeC3gTTQDdwmpTxbjbFoNJq5wxWrqlt6YbpRLRPQN6SUG6SUFwHbgL+o0jg0Go1mzlIVASClHLb9GgBkNcah0Wg0c5mq+QCEEF8BPggMAdeMctwdwB0AS5cunZrBaTQazRxASFmZzbcQ4kmgvcBbd0spf2o77rOAT0r5+bHOuWXLFrljx44yjlKj0WhmP0KInVLKLbmvV0wDkFJeV+KhDwCPAmMKAI1Go9GUj6r4AIQQq2y/vhPYX41xaDQazVymWj6Arwoh1mCEgZ4APlalcWg0Gs2cpSoCQEr5O9W4rkaj0WgyVMwJXAmEED0YGsNEaAV6yzicaqLnMv2YLfMAPZfpymTm0iGlbMt9cUYJgMkghNhRyAs+E9FzmX7MlnmAnst0pRJz0cXgNBqNZo6iBYBGo9HMUeaSAPhutQdQRvRcph+zZR6g5zJdKftc5owPQKPRaDTZzCUNQKPRaDQ2tACYoQghRLXHUC70XKYfs2UemtGZEwJACOEy/58VN7UQwg/U2X6fsfPSc5l+zJZ5aMZmVgsAIcRbhBA/Br4ihFgmZ7DDQxgEhBDfBF4FviGE+B2AmTav2TQXAHMuf8sMnsts+5sUQgjxW0KIJdUex2Qp5zxmrQAQQswHPofRcSwFfEkIcW11RzUxhBAu8yFcC5wPXA78GPgTIcSV5jEzYpcmhGg053I+M38uHUKIGmAlM3guQog62/21nhk6j2IIIdYJIV4C/i9wnxDiI9Ue00SoxDxmrQAANgFuKeX3gS8DvwHeJYRoru6wSkcIcaUQ4gHgc0KIDmAz8IKUskdK+UvgR8BX1OHVGmepCCHuBX4qhKgFLgBenolzEUJsEEL8FKOM+RrgSmbg30UIcZUQ4iHgfiHEdcBbgWdn2jxK4Cpgp5TycuCrwG8LIS6v8pgmQtnnMZsFwG4gZpp+IsBLQBK4vrrDGhshhEsI8Q3g74DHgRbgT4Fa4G3qOCnlt4DVQoiNUsr0dN2lCSHUfeYBnMANwBFsf4uZMheT24EdUsoLpJSvAl0YcwJmxlyEEJdhbIx+AjwGvAdwA+9Qx8yEeeQihFgkhPiGEOKTQoiF5sshICmEcEopnwLeAK42rQTTEiHEYiHEF4QQ/0sIsdh8eYQyz2M2C4Aw8BqZB/MwcAjoEEI4qzaqEpBSJoFfAL8tpfw3DJVvA/AQsFAIcbXt8PuAD5ufm5a2WnPx6MDo//xV4APAc0CLEMLeDnTaz0UIcR6wUEr5ZfP3jVLKHwGLcnZj030uq4CIlPI/gR8ATVLKv8KYx5W246b7PCyEEEsxtDInsAj4ByFEOxABhoEO89BtGPMv1LGwaph+GK8wuiQ+hfG8XAR8TQjRhGHKHqKM85jNAmAIeAXYKoRok1IGgXqMGz01A3Yzz0spTwshfFLKQSCBcWN/G0MgKI4Cp2fAfPowKhkeBwYxTHTfBP7YdsxMmMtRYKMQ4r2mGegfhRB/CpwC/ijnuOk8l6eANUKIb2GMdbkQ4nbgBeAztuOm9TyE0VdEsQx4Wkr5x1LKzwKngS8AT2AIhPMBpJT/A8wHlprnqPo6KIRYKw1iGJvX66SU/xv4OobwWgvsxFjsyzaPqk+8Upi7lV9gfJlfMF+uwVAHp/1uRkoZVf+bu84moMtUy5NCiL8QQtwIfAQ4Nt3nA1wHnJFSvgHswWgFehkQmElzMcd2P4ZJ7m+AtwDdQD+wVAjx50KIm5jmc5FSnsEw98wH3gu8HViMMY82IcT/mc7zEEKsF0I8A7yhopUwdsYdtsM+D/weIIFdwJuFEJvM9/ZhLpxSyvTUjDof2zxet83jPzGEFxibprVAt5TyEHAMuLRc85i1AgBAStmHsfg3CSFeAa7GWHhmGv8f8KQSCsCfAWeA/w18X0r546qNrHR6gAuFEI8An8BwKv4Qw6cxgLGgzpS5/Bq4EOOhTAIvYuyUv46xgP4JM2MuJzEEwA4p5QDwMhDD0DCn+zxiGCaQPyDTUfBB4HJzw4SpOf8cuBP4fxhWgXuEEP+EIfwen+IxF8I+jzsApJTdNoHbah6TMn//VwwNuizzmBO1gIQQbqBNSnm22mMZD6azJyWE+CrwOob28hHgK1LKF6s7uvFh7li+gmF6uBd4N/BOKeW7qzqwCSKE+A6QlFJ+UgjxaWCrlPK91R7XeBBCtGA4gg9KKf/OtD2vklLeXuWhlYQQoh5jI/EL4PNSyl+Yf5c6KeX7zWPeidF3/A7TF/VuYAXwb1LK7mqN3U7OPP6PlPJXQgivlDImhPgg8G4p5c05n/kdYDmTnMecEAAzGfPmOI1hZtgFfFdK+WR1RzV+hBDCbkYQQqwAaqSUe3LfmwkIIRoxtJYrgSDwJSnlSzNpLqZd/+3A3RiZv+eAu6WUr8ywefwxcLWU8p1CiHnAs8CdUsrHhBBfwNDU/rGqgywBcx5XSinfJYRwmALr6xgawjDwh8D9Usqny3bNGfI3nrMIIeowTD4/NkMOZzTCSGpLVnsc5UIIMW+67CQnihBiOeCQUh6p9lgmghBiAfDfwEellDuFEO8H3oRh8k0Bn5ZSPl/FIZaEbR63SSlfMzcZL2KET58D7pVS/qCs19QCQKPRzFSUpiKEuBMjZPIBDKfvs8CbpJTPVXN8pZIzjw0YgQZLgA8Bfy2l/EUlruuqxEk1Go1mKjAXTYERwfQB4DzgY2bAxIxY/CFvHh/CiPz5gDQqGVQMLQA0Gs1M5yZgIbBZSvlKtQczCaZ8HtoEpNFoZjQzyWE9GtWYhxYAGo1GM0eZ1YlgGo1GoymOFgAajUYzR9ECQKPRaOYoWgBoNBrNHEULAI1Go5mjaAGg0Wg0c5T/H3X9gnlGTUVCAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import seaborn as sns\n",
    "\n",
    "anomalies = test_score_df[test_score_df.anomaly == True]\n",
    "\n",
    "plt.plot(\n",
    "  range(test_dataset.data_len), \n",
    "  test_score_df['t'], \n",
    "  label='value'\n",
    ");\n",
    "\n",
    "sns.scatterplot(\n",
    "  anomalies.index,\n",
    "  anomalies.t,\n",
    "  color=sns.color_palette()[3],\n",
    "  s=52,\n",
    "  label='anomaly'\n",
    ")\n",
    "\n",
    "plt.plot(\n",
    "  range(len(test_score_df['y'])),\n",
    "  test_score_df['y'],\n",
    "  label='y'\n",
    ")\n",
    "\n",
    "plt.xticks(rotation=25)\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate the window-based anomalies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "start_end = []\n",
    "state = 0\n",
    "for idx in test_score_df.index:\n",
    "    if state==0 and test_score_df.loc[idx, 'y']==1:\n",
    "        state=1\n",
    "        start = idx\n",
    "    if state==1 and test_score_df.loc[idx, 'y']==0:\n",
    "        state = 0\n",
    "        end = idx\n",
    "        start_end.append((start, end))\n",
    "\n",
    "for s_e in start_end:\n",
    "    if sum(test_score_df[s_e[0]:s_e[1]+1]['anomaly'])>0:\n",
    "        for i in range(s_e[0], s_e[1]+1):\n",
    "            test_score_df.loc[i, 'anomaly'] = 1\n",
    "            \n",
    "actual = np.array(test_score_df['y'])\n",
    "predicted = np.array([int(a) for a in test_score_df['anomaly']])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate measurement scores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True Positive\t 7\n",
      "True Negative\t 105\n",
      "False Positive\t 1\n",
      "False Negative\t 8\n",
      "Accuracy\t 0.9256198347107438\n",
      "Precision\t 0.875\n",
      "Recall\t 0.4666666666666667\n",
      "f-measure\t 0.608695652173913\n",
      "cohen_kappa_score\t 0.5717656311443178\n",
      "auc\t 0.728616352201258\n",
      "roc_auc\t 0.728616352201258\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from sklearn.metrics import cohen_kappa_score\n",
    "from sklearn.metrics import roc_curve, auc, roc_auc_score\n",
    "\n",
    "predicted = np.array(predicted)\n",
    "actual = np.array(actual)\n",
    "\n",
    "tp = np.count_nonzero(predicted * actual)\n",
    "tn = np.count_nonzero((predicted - 1) * (actual - 1))\n",
    "fp = np.count_nonzero(predicted * (actual - 1))\n",
    "fn = np.count_nonzero((predicted - 1) * actual)\n",
    "\n",
    "print('True Positive\\t', tp)\n",
    "print('True Negative\\t', tn)\n",
    "print('False Positive\\t', fp)\n",
    "print('False Negative\\t', fn)\n",
    "\n",
    "accuracy = (tp + tn) / (tp + fp + fn + tn)\n",
    "precision = tp / (tp + fp)\n",
    "recall = tp / (tp + fn)\n",
    "fmeasure = (2 * precision * recall) / (precision + recall)\n",
    "cohen_kappa_score = cohen_kappa_score(predicted, actual)\n",
    "false_positive_rate, true_positive_rate, thresholds = roc_curve(actual, predicted)\n",
    "auc_val = auc(false_positive_rate, true_positive_rate)\n",
    "roc_auc_val = roc_auc_score(actual, predicted)\n",
    "\n",
    "print('Accuracy\\t', accuracy)\n",
    "print('Precision\\t', precision)\n",
    "print('Recall\\t', recall)\n",
    "print('f-measure\\t', fmeasure)\n",
    "print('cohen_kappa_score\\t', cohen_kappa_score)\n",
    "print('auc\\t', auc_val)\n",
    "print('roc_auc\\t', roc_auc_val)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import winsound\n",
    "# frequency = 200  # Set Frequency To 2500 Hertz\n",
    "# duration = 2000  # Set Duration To 1000 ms == 1 second\n",
    "# winsound.Beep(frequency, duration)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
